基于EDEM的青稞接觸參數(shù)仿真標(biāo)定*

丁文波，朱繼平，陳偉，張晉，袁棟，丁艷

(農(nóng)業(yè)農(nóng)村部南京農(nóng)業(yè)機(jī)械化研究所，南京市，210014)

0 引言

青稞作為西藏等地區(qū)人民的主糧，在高海拔、寒冷地區(qū)的糧食作物種植中，具有不可替代地位。青稞生產(chǎn)機(jī)械化程度低，已經(jīng)影響到青稞產(chǎn)業(yè)的發(fā)展和當(dāng)?shù)厝嗣竦纳钯|(zhì)量，提高青稞生產(chǎn)的機(jī)械化水平具有重要的現(xiàn)實(shí)意義[1-2]。

青稞機(jī)械化播種作為其機(jī)械化生產(chǎn)的重要環(huán)節(jié)，受限于排種器的研究發(fā)展，排種器的研究難點(diǎn)主要在于青稞種子在排種器中的運(yùn)動(dòng)及受力情況復(fù)雜，難以進(jìn)行深入研究。而EDEM離散元仿真技術(shù)已經(jīng)廣泛應(yīng)用于農(nóng)業(yè)領(lǐng)域[3-4]，能夠有效直觀的模擬機(jī)構(gòu)的工作原理，有利于高效便捷的對(duì)各種機(jī)構(gòu)進(jìn)行深入的研究[5-8]。進(jìn)行EDEM離散元仿真研究前需要建立物料的三維模型，設(shè)置物料的本征參數(shù)和接觸參數(shù)。建立精確的顆粒三維模型和準(zhǔn)確的設(shè)置仿真參數(shù)，可以提高仿真的精度。獲取這些參數(shù)有直接測量法[9-10]和虛擬標(biāo)定[11]兩種方法，通常一些難以直接測量的材料接觸參數(shù)需要通過虛擬標(biāo)定的方法進(jìn)行標(biāo)定。

目前對(duì)主要糧食作物如水稻、小麥、大豆、玉米等都已進(jìn)行相關(guān)的研究。張榮芳等[12]利用EDEM研究了不同填充球半徑的水稻種子模型對(duì)顆粒間的動(dòng)力學(xué)響應(yīng)特性的影響，尋找種子模型最佳的填充球顆粒數(shù)量，結(jié)果表明填充顆粒半徑為0.21 mm，仿真時(shí)長與仿真精度最優(yōu)；張濤等[13]對(duì)大豆排種過程的離散元仿真參數(shù)進(jìn)行研究，標(biāo)定了大豆種子的本征參數(shù)、大豆種子與排種器的接觸參數(shù)；崔濤等[14]通過離散元分析軟件EDEM仿真試驗(yàn)對(duì)玉米種子休止角試驗(yàn)所得試驗(yàn)數(shù)據(jù)進(jìn)行了驗(yàn)證，為排種器的設(shè)計(jì)與仿真分析提供理論依據(jù)。

目前對(duì)于青稞這方面的研究很少，本文以青稞藏青2000為研究對(duì)象，根據(jù)種子的形狀特征采用分段取截面的方法建立其三維模型，并通過臺(tái)架試驗(yàn)測定其部分參數(shù)，為其在仿真試驗(yàn)中的取值提供依據(jù)，最后以青稞休止角為指標(biāo)，通過響應(yīng)面法標(biāo)定青稞的主要接觸參數(shù)，以期為青稞后續(xù)的播種研究提供參考。

1 青稞的基本參數(shù)與三維模型

1.1 青稞的基本參數(shù)

青稞藏青2000由西藏自治區(qū)農(nóng)牧科學(xué)院農(nóng)業(yè)研究所于1995年以(藏青320×拉薩白青稞)與(喜拉19×昆侖164)復(fù)合雜交選育而成的青稞新品系，穗長方形、四棱長芒，黃白粒，具有良好的增產(chǎn)效果[15-16]。將經(jīng)過篩選的青稞種子隨機(jī)抽取五組測得平均千粒重46.5 g，含水率9.53%，密度10.2 g/cm3。隨機(jī)選取100粒種子使用游標(biāo)卡尺測得的平均長度7.75 mm、寬度3.76 mm、厚度2.67 mm。

1.2 青稞種子三維模型的建立

種子的三維模型與真實(shí)形狀越接近，仿真的結(jié)果就會(huì)越準(zhǔn)確，目前對(duì)于不規(guī)則的種子主要采取激光掃描的方式獲取精確的三維模型，如于慶旭等[17]利用逆向工程技術(shù)得到三七種子的精確三維輪廓模型；劉彩玲等[18]基于三維激光掃描獲取水稻種子精確的三維模型。對(duì)于小麥、青稞等近似橢球形，具有對(duì)稱性的種子主要采用近似法建模[19]。精確的三維建?？商岣叻抡娴木?，但精確建模難度大成本高。青稞藏青2000種子如圖1(a)所示，與小麥外形相似，故采用近似法建模。為了提高建模的精度，隨機(jī)抽取50粒種子，在長度方向進(jìn)行等分獲得等分截面如圖1(b)所示，測量其各個(gè)截面的寬度k、厚度h的尺寸參數(shù)，計(jì)算得到各截面寬、厚尺寸平均值如表1所示。再利用CATIA軟件中多截面實(shí)體功能，連接截面生成種子的一側(cè)實(shí)體，由種子的對(duì)稱性對(duì)一半實(shí)體進(jìn)行鏡像得到種子的三維模型如圖1(c)。

表1 青稞截面尺寸Tab. 1 Section size of highland barley

(a) 青稞種子 (b) 等分截面圖

2 臺(tái)架試驗(yàn)與仿真標(biāo)定

為了提高離散元仿真試驗(yàn)的可靠性，本研究利用臺(tái)架試驗(yàn)測量青稞種子與試驗(yàn)材料的靜摩擦因數(shù)以及青稞種子與試驗(yàn)材料的碰撞恢復(fù)系數(shù)，為這兩種離散元仿真參數(shù)的設(shè)置提供依據(jù)。對(duì)于不便直接測量的參數(shù)，采用臺(tái)架試驗(yàn)與仿真試驗(yàn)結(jié)合的方式進(jìn)行仿真標(biāo)定，以種子的休止角為檢驗(yàn)指標(biāo)進(jìn)行仿真試驗(yàn)，休止角采用圓筒提升法獲得[20]。試驗(yàn)材料選擇3D打印材料pla(聚乳酸)塑料[21]，方便為后續(xù)3D打印排種器進(jìn)行排種機(jī)理研究做準(zhǔn)備。

2.1 臺(tái)架試驗(yàn)

2.1.1 斜面法確定靜摩擦系數(shù)取值

通過斜面法測量青稞種子與pla材料的靜摩擦系數(shù)，試驗(yàn)材料有pla3D打印平板如圖2(a)所示(長150 mm、寬50 mm、厚2 mm)、青稞藏青2000種子、萬能角度尺、多功能鐵架臺(tái)等。為了防止青稞種子在試驗(yàn)中發(fā)生滾動(dòng)影響試驗(yàn)的結(jié)果，將2粒青稞種子粘在一起為一組，共設(shè)置五組如圖2(b)所示。試驗(yàn)的原理：將pla平板置于多功能鐵架臺(tái)的可轉(zhuǎn)動(dòng)鐵夾上，轉(zhuǎn)動(dòng)鐵夾使得平板處于水平位置，將粘好的種子放置上面，然后緩慢轉(zhuǎn)動(dòng)鐵夾使得平板逐漸傾斜，直到種子發(fā)生滑動(dòng)時(shí)停止轉(zhuǎn)動(dòng)鐵夾。用萬能角度尺測量此時(shí)的平板的傾斜角，每組重復(fù)5次，完成平板傾角的測定，計(jì)算平均值求得平板傾斜角α=27.04°。根據(jù)公式求得青稞與pla平板的最大靜摩擦系數(shù)x8=0.51。由于試驗(yàn)本身存在誤差，故x8取值應(yīng)該在0.51附近，取值范圍定為0.3～0.7具有可信度。

(a) pla 3D打印平板 (b) 青稞種子粘結(jié)處理

2.1.2 臺(tái)架試驗(yàn)確定恢復(fù)系數(shù)取值

對(duì)青稞—pla平板之間的碰撞恢復(fù)系數(shù)的測量是通過使青稞種子從一定的高度垂直落下，與放置于正下方呈45°放置的pla平板發(fā)生碰撞后，作拋物線運(yùn)動(dòng)，最后落入收集裝置[22]，恢復(fù)系數(shù)測定裝置如圖3(a)，通過測量種子初始位置N點(diǎn)與碰撞點(diǎn)O的垂直高度h1、碰撞點(diǎn)O與落點(diǎn)P的垂直高度h2、碰撞點(diǎn)與落點(diǎn)的水平距離s1等位移，代入運(yùn)動(dòng)學(xué)方程[23]求解得到恢復(fù)系數(shù)。重復(fù)試驗(yàn)10次記錄數(shù)據(jù)，通過計(jì)算求得碰撞恢復(fù)系數(shù)平均值x7=0.49。試驗(yàn)裝置原理圖如圖3(b)。

(a) 恢復(fù)系數(shù)測定裝置

2.1.3 圓筒提升試驗(yàn)確定休止角

采用圓筒提升法測量青稞休止角θ[24]。試驗(yàn)的材料有內(nèi)徑高度分別為21 mm、200 mm的玻璃圓筒、多功能鐵架臺(tái)、直徑厚度分別為80 mm和3 mm的pla圓板、萬能角度尺、游標(biāo)卡尺等。將pla圓板懸空置于鐵架上，然后將充滿青稞顆粒的圓筒口對(duì)準(zhǔn)圓板的中心倒置，再以0.02 m/s的速度緩慢提升圓筒，最后形成穩(wěn)定的顆粒堆。通過萬能角度尺測量得到顆粒堆休止角θ。重復(fù)10次試驗(yàn)求得平均休止角θ=25.1°，試驗(yàn)裝置如圖4。

圖4 圓筒提升試驗(yàn)裝置圖

2.2 仿真模型的建立及參數(shù)設(shè)置

2.2.1 青稞的離散元仿真模型

將種子的三維模型(圖1(c))導(dǎo)入EDEM軟件中，采用多球聚合模型[25]對(duì)種子模型進(jìn)行填充，通過若干個(gè)直徑不等的球形顆粒重疊堆積完成，得到的種子離散元模型如圖5所示，從圖中可以看出顆粒模型與種子實(shí)際外輪廓較為吻合。

圖5 青稞的離散元仿真顆粒模型

2.2.2 仿真參數(shù)設(shè)置

離散元顆粒堆積仿真試驗(yàn)所需要的基本參數(shù)有青稞和pla材料的本征參數(shù)及接觸參數(shù)。其中已知的本征參數(shù)如表2所示：青稞的密度由試驗(yàn)測得為1.02 g/cm3；通過查閱文獻(xiàn)[26]得到pla材料的泊松比0.29、剪切模量220 MPa、密度1.11 g/cm3。

表2 仿真參數(shù)的取值Tab. 2 Value of simulation parameters

未確定的仿真參數(shù)通過查找文獻(xiàn)[27-28]及臺(tái)架試驗(yàn)確定其取值的范圍如表2所示。通過臺(tái)架試驗(yàn)測得x8=0.51、x7=0.49。由于試驗(yàn)本身存在一定的誤差，故以試驗(yàn)測定的值為中心值擴(kuò)展，得到這兩個(gè)仿真參數(shù)的取值區(qū)間，具有可信度。仿真顆粒的分布按照標(biāo)準(zhǔn)正態(tài)分布設(shè)置，平均值Mean：1；標(biāo)準(zhǔn)差StdDev：0.05；由青稞種子顆粒的運(yùn)動(dòng)特點(diǎn)，接觸模型選擇“Hertz-Mindlin(noslip)”；仿真中FixedTimeStep采用20%的瑞利時(shí)步；仿真中網(wǎng)格尺寸CellSize取3倍最小球形單元尺寸。

2.2.3 圓筒提升仿真試驗(yàn)

在圓筒提升試驗(yàn)中影響試驗(yàn)指標(biāo)的因素有：青稞的泊松比x1、青稞的剪切模量x2、青稞—青稞的碰撞恢復(fù)系數(shù)x3、青稞—青稞的靜摩擦系數(shù)x4、青稞—青稞的滾動(dòng)摩擦系數(shù)x5、青稞—pla平板的滾動(dòng)摩擦系數(shù)x6、青稞—pla平板的碰撞恢復(fù)系數(shù)x7、青稞—pla平板的靜摩擦因數(shù)x8。試驗(yàn)的指標(biāo)：臺(tái)架試驗(yàn)中青稞的休止角θ。

利用三維軟件建立出圓筒提升試驗(yàn)中圓筒和pla圓板的三維模型，將其保存為Stp格式導(dǎo)入EDEM軟件中。在圓筒內(nèi)設(shè)置顆粒工廠生成顆粒，經(jīng)過預(yù)仿真試驗(yàn)確定以1 500個(gè)/s的速率生成顆粒，生成時(shí)間設(shè)置為1 s，可產(chǎn)生足夠的試驗(yàn)顆粒，顆粒生成后自由落下填充圓筒。在顆粒生成結(jié)束后，給圓筒添加一個(gè)向上的直線運(yùn)動(dòng)，速度為0.02 m/s，隨著圓筒的提升，顆粒緩慢從桶內(nèi)流出，在pla圓板上形成穩(wěn)定的顆粒堆，多余的顆粒則落入下方的容器中。仿真完成以后，進(jìn)入EDEM后處理界面，利用軟件自帶量角功能直接測量顆粒堆的休止角θ1，如圖6所示。

圖6 青稞顆粒休止角測量

2.3 仿真模型的響應(yīng)面試驗(yàn)設(shè)計(jì)

2.3.1 Plackett-Burman試驗(yàn)

通過Plackett-Burman試驗(yàn)篩選出顯著性試驗(yàn)因素，并非所有試驗(yàn)因素都對(duì)試驗(yàn)指標(biāo)顆粒堆積角有顯著的影響，對(duì)于影響不顯著的試驗(yàn)因素，無法進(jìn)行仿真分析。故應(yīng)用Design-Expert軟件進(jìn)行Plackett-Burman試驗(yàn)，仿真試驗(yàn)因素及水平如表3所示，篩選出x1～x8中對(duì)試驗(yàn)指標(biāo)有顯著影響的因素。仿真試驗(yàn)設(shè)計(jì)及結(jié)果如表4所示。

表3 Plackett-Burman試驗(yàn)因素及水平Tab. 3 Factors and levels of Plackett-Burman test

表4 Plackett-Burman 試驗(yàn)設(shè)計(jì)及結(jié)果Tab. 4 Design and results of Plackett-Burman test

通過Design-Expert的分析功能得到表5，從表5可以看出對(duì)青稞休止角影響最顯著的參數(shù)分別是x4、x5、x8是影響休止角大小的主要因素，其他的影響因素對(duì)休止角的影響不顯著。

表5 Plackett-Burman試驗(yàn)參數(shù)顯著性分析Tab. 5 Analysis of significance of parameters in Plackett-Burman test

2.3.2 最陡爬坡試驗(yàn)

采用最陡坡試驗(yàn)可以較為快速的確定顯著性因素的最佳取值范圍，為后續(xù)的響應(yīng)面試驗(yàn)提供取值依據(jù)，從而確定參數(shù)的最優(yōu)值。由前文中Plackett-Burman試驗(yàn)可知顯著性影響因素有x4、x5、x8，取值范圍分別為0.1～0.6、0～0.1、0.3～0.7，等分取六個(gè)梯度值較為合適。其他因素對(duì)試驗(yàn)影響不顯著，則都取中間值，分別取x1(0.125)、x2(5.5 MPa)、x3(0.5)、x6(0.05)、x7(0.5)。

最陡爬坡試驗(yàn)設(shè)計(jì)及結(jié)果如表6所示，結(jié)果表明隨著x4、x5、x8的增大仿真試驗(yàn)的顆粒休止角也不斷增大，仿真試驗(yàn)的顆粒休止角與臺(tái)架試驗(yàn)測得的休止角之間的相對(duì)誤差先減小后增大。在第二梯度時(shí)相對(duì)誤差最小，故最優(yōu)值在2號(hào)附近，則最優(yōu)的取值區(qū)間在1號(hào)和3號(hào)之間。故表7中1、2、3號(hào)取值可分別作為后面Box-Behnken試驗(yàn)中-1、0、1水平的取值。

表6 最陡爬坡試驗(yàn)及結(jié)果Tab. 6 Design and results of steepest ascent test

2.3.3 Box-Behnken試驗(yàn)及回歸模型

根據(jù)Plackett-Burman試驗(yàn)與最陡爬坡試驗(yàn)，確定了顯著性影響因素有x4、x5、x8及其在Box-Behnken試驗(yàn)中-1、0、1水平的取值。其他因素的取值與最陡爬坡試驗(yàn)中相同。采用Design-Expert軟件進(jìn)行Box-Behnken試驗(yàn)設(shè)計(jì)并建立回歸模型，共進(jìn)行17次仿真試驗(yàn)。Box-Behnken試驗(yàn)設(shè)計(jì)及結(jié)果如表7所示。根據(jù)模型得到其二次多項(xiàng)式方程如式(1)。

θ1=26.67+5.12x4+1.51x5-0.29x8-

0.63x4x5+0.15x4x8+0.15x5x8+

0.61x42-0.81x52+0.20x82

(1)

表7 Box-Behnken試驗(yàn)設(shè)計(jì)及結(jié)果Tab. 7 Design and results of Box-Behnken test

仿真試驗(yàn)的方差分析如表8所示，模型的p值小于0.01該回歸模型的擬合性好；x4、x5等因素的p值均小于0.01，說明這兩個(gè)因素都對(duì)顆粒休止角具有極顯著性的影響；x8的p值小于0.05說明對(duì)休止角具有顯著性影響；交互項(xiàng)x4x5、x42、x52的p值都小于0.01說明其對(duì)休止角也具有極顯著性影響，交互項(xiàng)x4x5對(duì)休止角的影響如圖7所示，x4、x5、都對(duì)休止角起正效應(yīng)作用；交互項(xiàng)x4x8、x5x8、x82的p值都大于0.05，對(duì)顆粒的休止角沒有顯著影響；決定系數(shù)R2=0.998 2，校正決定系數(shù)R2adj=0.996 0，二者都接近1，表明擬合方程具有較高的可靠性。

在保證模型顯著、失擬項(xiàng)不顯著的基礎(chǔ)上，剔除模型中的不顯著回歸項(xiàng)，對(duì)模型進(jìn)行優(yōu)化[29]得到新的回歸方程如式(2)。

θ1=26.76+5.12x4+1.51x5-0.29x8-

0.63x4x5+0.62x42-0.80x52

(2)

最后，將臺(tái)架試驗(yàn)測定的青稞休止角代入到Design-Expert中，求得3個(gè)顯著性參數(shù)的最優(yōu)值，x4(0.19)、x5(0.01)、x8(0.43)。

表8 Box-Behnken 試驗(yàn)設(shè)計(jì)二次多項(xiàng)式模型方差分析Tab. 8 ANOVA of quadratic polynomial model of Box-Behnken test

圖7 x4與x5交互作用

2.4 試驗(yàn)驗(yàn)證

將Box-Behnken試驗(yàn)得到的三個(gè)顯著性參數(shù)的最優(yōu)值代入到EDEM仿真模型中，進(jìn)行仿真驗(yàn)證試驗(yàn)。仿真參數(shù)設(shè)置：青稞—青稞的靜摩擦系數(shù)x4=0.19、青稞—青稞的滾動(dòng)摩擦系數(shù)x5=0.01、青稞—pla平板的靜摩擦因數(shù)x8=0.43，其他非顯著性因素與之前相同，重復(fù)5次重復(fù)試驗(yàn)得到顆粒的平均休止角θ1=25.26°，與實(shí)際青稞休止角25.1°相對(duì)誤差0.64%，相對(duì)誤差很小，說明標(biāo)定的參數(shù)可信，可以用于青稞離散元仿真，為后續(xù)青稞播種等研究提供參考。

3 結(jié)論

1) 通過Plackett-Burman試驗(yàn)，從8個(gè)試驗(yàn)因素中找出了主要影響試驗(yàn)指標(biāo)青稞休止角的三個(gè)顯著性因素：青稞—青稞的靜摩擦系數(shù)x4、青稞—青稞的滾動(dòng)摩擦系數(shù)x5、青稞—pla 3D打印材料的靜摩擦因數(shù)x8。其他因素為非顯著性因素，對(duì)試驗(yàn)指標(biāo)的影響不顯著。

2) 根據(jù)最陡爬坡試驗(yàn)確定3個(gè)顯著性參數(shù)的最優(yōu)值取值范圍，為后面進(jìn)行Box-Behnken試驗(yàn)提供高低水平的取值，以此為基礎(chǔ)建立顯著性參數(shù)與休止角之間的二次回歸模型并對(duì)其進(jìn)行優(yōu)化，試驗(yàn)表明：3個(gè)顯著性因素對(duì)顆粒休止角具有顯著的影響，此外青稞與青稞的靜摩擦系數(shù)的二次項(xiàng)、青稞與青稞的滾動(dòng)摩擦系數(shù)的二次項(xiàng)、二者的交互項(xiàng)對(duì)顆粒休止角也有顯著影響。

3) 取響應(yīng)面試驗(yàn)得到青稞—青稞的靜摩擦系數(shù)x4為0.19、青稞—青稞的滾動(dòng)摩擦系數(shù)x5為0.01、青稞—pla平板的靜摩擦因數(shù)x8為0.43。其余非顯著性參數(shù)取值：青稞的泊松比x1為0.18、青稞的剪切模量x2為5.5 MPa、青稞—青稞的碰撞恢復(fù)系數(shù)x3為0.5、青稞—pla平板的滾動(dòng)摩擦系數(shù)x6為0.05、青稞—pla平板的撞恢復(fù)系數(shù)x7為0.5。將以上參數(shù)導(dǎo)入EDEM中進(jìn)行仿真驗(yàn)證試驗(yàn)與實(shí)際臺(tái)架試驗(yàn)的結(jié)果進(jìn)行比較，平均休止角θ1=25.26°與實(shí)際青稞休止角25.1°相對(duì)誤差0.64%，誤差很小，表明本研究標(biāo)定的青稞離散元仿真參數(shù)具有較高的可靠性，可以為后續(xù)的青稞仿真研究提供參考。