基于EDEM的青稞接觸參數仿真標定*

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

(農業農村部南京農業機械化研究所，南京市，210014)

0 引言

青稞作為西藏等地區人民的主糧，在高海拔、寒冷地區的糧食作物種植中，具有不可替代地位。青稞生產機械化程度低，已經影響到青稞產業的發展和當地人民的生活質量，提高青稞生產的機械化水平具有重要的現實意義[1-2]。

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

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

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

1 青稞的基本參數與三維模型

1.1 青稞的基本參數

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

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

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

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

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

2 臺架試驗與仿真標定

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

2.1 臺架試驗

2.1.1 斜面法確定靜摩擦系數取值

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

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

2.1.2 臺架試驗確定恢復系數取值

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

(a) 恢復系數測定裝置

2.1.3 圓筒提升試驗確定休止角

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

圖4 圓筒提升試驗裝置圖

2.2 仿真模型的建立及參數設置

2.2.1 青稞的離散元仿真模型

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

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

2.2.2 仿真參數設置

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

表2 仿真參數的取值Tab. 2 Value of simulation parameters

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

2.2.3 圓筒提升仿真試驗

在圓筒提升試驗中影響試驗指標的因素有：青稞的泊松比x1、青稞的剪切模量x2、青稞—青稞的碰撞恢復系數x3、青稞—青稞的靜摩擦系數x4、青稞—青稞的滾動摩擦系數x5、青稞—pla平板的滾動摩擦系數x6、青稞—pla平板的碰撞恢復系數x7、青稞—pla平板的靜摩擦因數x8。試驗的指標：臺架試驗中青稞的休止角θ。

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

圖6 青稞顆粒休止角測量

2.3 仿真模型的響應面試驗設計

2.3.1 Plackett-Burman試驗

通過Plackett-Burman試驗篩選出顯著性試驗因素，并非所有試驗因素都對試驗指標顆粒堆積角有顯著的影響，對于影響不顯著的試驗因素，無法進行仿真分析。故應用Design-Expert軟件進行Plackett-Burman試驗，仿真試驗因素及水平如表3所示，篩選出x1～x8中對試驗指標有顯著影響的因素。仿真試驗設計及結果如表4所示。

表3 Plackett-Burman試驗因素及水平Tab. 3 Factors and levels of Plackett-Burman test

表4 Plackett-Burman 試驗設計及結果Tab. 4 Design and results of Plackett-Burman test

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

表5 Plackett-Burman試驗參數顯著性分析Tab. 5 Analysis of significance of parameters in Plackett-Burman test

2.3.2 最陡爬坡試驗

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

最陡爬坡試驗設計及結果如表6所示，結果表明隨著x4、x5、x8的增大仿真試驗的顆粒休止角也不斷增大，仿真試驗的顆粒休止角與臺架試驗測得的休止角之間的相對誤差先減小后增大。在第二梯度時相對誤差最小，故最優值在2號附近，則最優的取值區間在1號和3號之間。故表7中1、2、3號取值可分別作為后面Box-Behnken試驗中-1、0、1水平的取值。

表6 最陡爬坡試驗及結果Tab. 6 Design and results of steepest ascent test

2.3.3 Box-Behnken試驗及回歸模型

根據Plackett-Burman試驗與最陡爬坡試驗，確定了顯著性影響因素有x4、x5、x8及其在Box-Behnken試驗中-1、0、1水平的取值。其他因素的取值與最陡爬坡試驗中相同。采用Design-Expert軟件進行Box-Behnken試驗設計并建立回歸模型，共進行17次仿真試驗。Box-Behnken試驗設計及結果如表7所示。根據模型得到其二次多項式方程如式(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試驗設計及結果Tab. 7 Design and results of Box-Behnken test

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

在保證模型顯著、失擬項不顯著的基礎上，剔除模型中的不顯著回歸項，對模型進行優化[29]得到新的回歸方程如式(2)。

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

0.63x4x5+0.62x42-0.80x52

(2)

最后，將臺架試驗測定的青稞休止角代入到Design-Expert中，求得3個顯著性參數的最優值，x4(0.19)、x5(0.01)、x8(0.43)。

表8 Box-Behnken 試驗設計二次多項式模型方差分析Tab. 8 ANOVA of quadratic polynomial model of Box-Behnken test

圖7 x4與x5交互作用

2.4 試驗驗證

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

3 結論

1) 通過Plackett-Burman試驗，從8個試驗因素中找出了主要影響試驗指標青稞休止角的三個顯著性因素：青稞—青稞的靜摩擦系數x4、青稞—青稞的滾動摩擦系數x5、青稞—pla 3D打印材料的靜摩擦因數x8。其他因素為非顯著性因素，對試驗指標的影響不顯著。

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

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