花生莢果離散元仿真參數(shù)標定*

徐效偉，魏海，顏建春，鮑國丞，杜元杰，謝煥雄

(1. 農(nóng)業(yè)農(nóng)村部南京農(nóng)業(yè)機械化研究所，南京市，210014; 2. 江蘇省農(nóng)業(yè)機械試驗鑒定站，南京市，210017)

0 引言

花生(peanut)，原名落花生，又名“長生果”“泥豆”“番豆”?；ㄉ且环N重要的油料作物，同時也是優(yōu)質(zhì)蛋白質(zhì)來源。2019年全國花生種植面積約為4 633 khm2，產(chǎn)量為1 752 kt，花生種植面積和總產(chǎn)量常年居于世界前列[1-8]。近年來花生生產(chǎn)機械化率逐年提高，但相關機械較復雜[9-12]，通過傳統(tǒng)試驗進行前期的試驗驗證，費時費工，基于離散元法的EDEM軟件在工農(nóng)業(yè)中應用越來越廣泛，因此運用EDEM軟件進行前期仿真試驗，可以有效縮短研發(fā)周期、降低研發(fā)成本。

目前，離散元仿真分析參數(shù)標定的研究對象多為肥料、土壤和種子等。韓樹杰等[13]采用物理試驗與仿真試驗相結(jié)合的方法，以堆積角為響應值對散體廄肥離散元參數(shù)進行標定。邢潔潔等[14]對海南地區(qū)磚紅土壤離散元接觸參數(shù)進行仿真試驗標定，以實際堆積角為目標值，進行仿真試驗，獲得最佳參數(shù)組合。于慶旭等[15]采用逆向工程技術對三七種子進行三維建模，并利用實際試驗與仿真試驗相結(jié)合的方法對三七種子的相關物理參數(shù)進行了標定，在后續(xù)的臺架試驗驗證中，試驗效果與標定后仿真結(jié)果一致。

花生莢果相關物理參數(shù)研究主要集中在摩擦系數(shù)和恢復系數(shù)等物理量。目前對于某一品種花生基本物理參數(shù)研究較少，離散元仿真分析在花生機械研發(fā)中應用較少。在花生莢果仿真物理模型參數(shù)標定及離散元仿真分析方面應用鮮有報道[16]。

本文以宛花2號花生莢果為研究對象，通過物理試驗和仿真試驗分析相結(jié)合的方法確定花生莢果基本物理參數(shù)和接觸力學參數(shù)，并建立花生莢果離散元仿真三維模型[17]。使用EDEM軟件(離散元法)對花生莢果顆粒堆積過程進行仿真試驗，運用MATLAB軟件進行圖像處理測定花生莢果的堆積角，利用篩選試驗、最陡爬坡試驗及響應面試驗確定花生莢果仿真參數(shù)間的最佳組合，以期為花生莢果采摘、烘干和脫殼等作業(yè)仿真提供參數(shù)依據(jù)。

1 材料和方法

1.1 試驗材料與實驗儀器

本文選用河南省主要種植品種宛花2號花生莢果為試驗物料，試驗所用花生莢果經(jīng)初步清理和篩選，無損傷，無霉變，該花生莢果采用兩段式收獲方式(先挖掘晾曬，再撿拾收獲)。試驗主要用實驗儀器有DWD型電子萬能試驗機(額定負荷5 kN、準確度等級0.5、位移分辨率0.01 mm、加載速度0.01～500 mm/min)、數(shù)顯游標卡尺(精度0.01 mm)、自制摩擦系數(shù)測量儀、單軸傾角盒(精度0.05°)、DGF30/7-IA型電熱鼓風干燥箱(溫度0 ℃～300 ℃，電壓220 V)、上皿電子天平(量程100 g，精度0.000 1 g)、FASTEC IMAGING生產(chǎn)的HiSpec5型高速攝像機、50 mL量筒(精度1 mL)、自制鐵盒(容積1.047 625×10-3m3)。

1.2 物理參數(shù)測定

1.2.1 基本物理參數(shù)測定

花生莢果屬于散粒體物料，其基本物理特性包括花生莢果外形尺寸(長X×寬Y×厚Z)、密度、含水率、容重、泊松比、彈性模量和剪切模量[18]。隨機選取200顆花生莢果平均分為5組。利用量筒、上皿電子天平和數(shù)顯游標卡尺測定花生莢果的外形尺寸和密度。根據(jù)GB/T 20264—2006《糧食、油料水分兩次烘干測定法》，利用烘干箱和電子天平進行含水率測量[19-20]。用自制鐵盒測定花生莢果的容重。

通過20次重復試驗，測得花生莢果的外形尺寸(長X×寬Y×厚Z)、密度、容重和含水率如表1所示。

表1 花生莢果基本物理參數(shù)Tab. 1 Basic physical parameters of peanut pods

1.2.2 泊松比

宛花2號花生莢果形狀不一，且外形不規(guī)則。隨機選取10顆花生莢果，記錄其長X和寬Y原始尺寸。利用DWD型電子萬能試驗機對整個花生莢果進行壓力變形試驗，通過測量花生莢果加載前后開裂極限處寬度方向和長度方向的變形量計算泊松比。試驗過程中采用平板壓頭以10 mm/min速度沿花生莢果長度方向加載，當花生莢果出現(xiàn)破裂時停止加載。由電子萬能試驗機記錄其長度方向變形量，數(shù)顯游標卡尺記錄花生莢果在寬度方向載荷開裂極限處的變形量。通過式(1)計算泊松比，結(jié)果取平均值。

(1)

式中：μ——花生莢果泊松比；

εX——寬度方向應變；

εY——長度方向應變；

ΔX——長度方向變形量，mm；

ΔY——寬度方向變形量，mm。

試驗重復10次取平均值，測得花生莢果的泊松比為0.4。

1.2.3 彈性模量和剪切模量計算

花生莢果不是規(guī)則的圓柱體，若按照常規(guī)壓縮試驗法壓縮整個花生莢果測定彈性模量誤差較大，試驗的難點在于確定壓縮過程中花生莢果橫截面面積?；ㄉv果彈性模量主要與莢果外殼相關，因此本試驗選取形狀近圓柱的花生莢果，并從中截取近圓柱體部分作為試樣并量取試樣高度，分別量取試樣上、中、下三處內(nèi)外直徑，計算面積后取平均值作為花生莢果壓縮過程中的橫截面積。本文采用管件壓縮試驗測定彈性模量，試驗時將試樣放置在電子萬能試驗機平板上，并設定加載速度和加載時間分別為12 mm/min和10 s，使用平板壓頭，沿試樣軸向壓縮10 s后停機，從軟件中獲取試樣壓縮過程中載荷、位移相關數(shù)據(jù)。試驗重復10次取平均值，通過式(2)計算得花生莢果試樣彈性模量為17.46±5.72 MPa。試樣壓縮如圖1所示。

(2)

式中：E——花生莢果彈性模量，MPa；

FN——壓縮過程中的最大力，N；

l——試樣原始長度，mm；

Δl——壓縮過程中最大力下的變形量，mm；

A——試樣壓縮受力橫截面積，mm2。

圖1 花生莢果試樣壓縮試驗

根據(jù)試驗測定的泊松比和彈性模量，通過式(3)計算可得花生莢果剪切模量為6.5 MPa。

(3)

式中：G——花生莢果剪切模量，MPa。

1.3 接觸參數(shù)測定

利用EDEM軟件進行離散元仿真試驗時所需要的接觸參數(shù)有花生莢果之間和花生莢果與鋼板之間的碰撞恢復系數(shù)、靜摩擦系數(shù)和滾動摩擦系數(shù)，通過實際物理試驗對相關基礎參數(shù)進行測定，為離散元仿真試驗提供參考。

1.3.1 碰撞恢復系數(shù)測定

碰撞恢復系數(shù)是衡量物體碰撞后恢復到原來形狀的能力[18]。碰撞恢復系數(shù)定義為碰撞前后兩物體沿接觸處法線方向上的分離速度與接近速度之比[21-22]。即

(4)

式中：e——物體的恢復系數(shù)；

v1、v2——碰撞后物體1和2的速度，m/s；

v10、v20——碰撞前物體1和2的速度，m/s。

(5)

式中：h——花生莢果彈起后最大高度，mm；

H——花生莢果自由落體的高度，mm。

圖2 恢復系數(shù)測試原理圖

根據(jù)上述原理，本文采用HiSpec5型高速攝像機采集花生莢果跌落碰撞視頻與照片，攝像機設置為150幀/s，試驗裝置主要包括高速攝像機、跌落架和間距為3 mm的方格背景紙，攝像機鏡頭與方格背景紙之間的距離為50 cm，跌落碰撞試驗裝置如圖3所示。本文分別以Q235鋼板和花生莢果整齊排列的莢果板作為碰撞接觸材料，進行花生莢果跌落碰撞試驗，用于測取花生莢果與Q235鋼板和花生莢果之間的碰撞恢復系數(shù)。花生莢果從距離碰撞接觸材料40 cm處自由下落，與被測接觸材料碰撞后彈起，利用高速攝像機捕捉花生莢果試驗過程。

圖3 跌落碰撞試驗裝置

利用高速攝像機控制軟件，將跌落碰撞試驗的視頻以圖片形式導出，截取跌落碰撞過程圖片，找出花生莢果碰撞后到達最高點處的圖片，并數(shù)出背景方格紙的方格數(shù)，由式(5)計算得出花生莢果的碰撞恢復系數(shù)[23]。

利用上述方法分別測取花生莢果與篩板、花生莢果之間的碰撞恢復系數(shù)，每個試驗重復10次，測取結(jié)果取平均值，花生莢果與Q235鋼板之間恢復系數(shù)0.275±0.064，花生莢果之間恢復系數(shù)0.287±0.043。

1.3.2 靜摩擦系數(shù)測定

本文采用靜力學原理對花生莢果摩擦系數(shù)進行測定，花生莢果在自制摩擦系數(shù)測量儀上的靜力學分析如圖4所示。由靜力學原理可知，在花生莢果受力滿足式(6)[24]；當斜面傾角逐漸增大到使花生莢果有下滑趨勢時，滿足式(7)。

圖4 摩擦系數(shù)測量原理示意

(6)

μ1=f/N=tanθ

(7)

式中：θ——花生莢果的靜摩擦角，(°)；

μ1——花生莢果靜摩擦系數(shù)；

m——花生莢果質(zhì)量，kg；

g——重力加速度，9.8 N/s2。

此時斜面傾角為花生莢果最大靜摩擦角，利用該方法可以測定花生莢果的靜摩擦系數(shù)。本試驗利用自制摩擦系數(shù)測量儀(圖5)和單軸傾角盒測量所需花生莢果與各接觸材料間的靜摩擦系數(shù)。測量花生莢果與Q235鋼板的靜摩擦系數(shù)時，將處理好的花生莢果平放在桌面上，用大頭針固定花生莢果，先將單軸傾角盒放置于摩擦系數(shù)測量儀的鉸接軸附近，避免影響花生滑動，再將花生莢果放置于摩擦系數(shù)測量儀上，保證測試花生莢果與接觸材料面充分接觸，緩慢轉(zhuǎn)動搖把增大斜面傾斜角度，當花生莢果沿斜面突然下滑時，搖把停止轉(zhuǎn)動，記錄此時單軸傾角盒的數(shù)值。該值即為待測花生莢果在測試板上的靜摩擦角。測定花生莢果間靜摩擦系數(shù)時，先將花生莢果粘貼在平面板上，制成莢果板，然后重復上述試驗過程，記錄靜摩擦角，計算莢果間的靜摩擦系數(shù)。每組試驗重復20次計算后取平均值，得到花生莢果與Q235鋼板間靜摩擦系數(shù)平均值為0.583±0.015，花生莢果間靜摩擦系數(shù)平均值為0.823±0.041。

圖5 摩擦系數(shù)斜面儀

1.3.3 動摩擦系數(shù)測定

動摩擦系數(shù)由自制摩擦系數(shù)測量儀測定，花生莢果在斜面抬升過程中出現(xiàn)滾動時斜面傾角為花生莢果動摩擦角[25]，每組試驗重復20次，計算后取平均值，花生莢果與Q235鋼板間動摩擦系數(shù)平均值為0.104±0.024，花生莢果間動摩擦系數(shù)平均值為0.283±0.043。由于花生莢果表面不規(guī)則，且表面粗糙，試驗過程中系數(shù)變化較大，因此在后續(xù)的仿真試驗過程中利用堆積角試驗對動摩擦系數(shù)進行標定。

1.3.4 堆積角測定

本試驗選用Q235鋼圓筒，根據(jù)花生莢果尺寸，確定其內(nèi)徑和高度分別為15 cm和45 cm。試驗時圓筒垂直放置在Q235鋼板上，花生莢果經(jīng)過初步篩選后填充滿整個圓筒。隨后將圓筒以0.05 m/s的速度勻速提升[26]，底板上形成近似錐型的花生莢果堆，試驗重復10次。每次試驗后拍攝莢果堆正視圖像，利用Matlab對花生莢果堆積圖像進行灰度處理、二值化處理、圖像邊界提??；利用Origin對邊界數(shù)據(jù)進行擬合，得到花生莢果堆積角平均值為31.63°，如圖6所示?；ㄉv果單側(cè)堆積角邊界擬合如圖7所示。

(a) 原圖像 (b) 灰度圖像 (c) 二值化圖像

圖7 單側(cè)堆積角邊界擬合

2 花生莢果離散元模型建立及仿真參數(shù)標定

2.1 離散元仿真模型建立

利用Inventor軟件按照花生莢果實際測量所得尺寸建立三維模型(顆粒模型)，將花生莢果模型轉(zhuǎn)成.stl格式后倒入離散元仿真軟件EDEM中。對于仿真模型的建立，其尺寸和形狀對仿真結(jié)果影響較小，顆?！w粒、顆粒—鋼板之間的相互作用對仿真結(jié)果影響較大[17]。本文通過EDEM軟件中單球形顆粒對花生莢果模型進行填充，建立花生莢果的仿真模型，如圖8所示。

圖8 花生莢果仿真三維模型

在EDEM仿真模擬試驗中采用Hertz-Mindin模型來模擬花生莢果在干燥倉內(nèi)的運動特性及物料之間的相互作用。花生莢果仿真模擬堆積角試驗模型按照實際堆積角測定裝置進行繪制，如圖9所示。在導管上方口徑處建立顆粒工廠，用于生成花生莢果顆粒，顆粒采用動態(tài)生成方式，生成速率為4 800粒/s，一共生成1 200粒，顆粒尺寸采用固定形式。仿真試驗總時間為5 s，時間步長為2.88×10-7s，網(wǎng)格尺寸大小為最小顆粒半徑的3倍[27]。

圖9 堆積角測定裝置

仿真開始時，花生莢果顆粒從導管上方口徑處的虛擬顆粒平面生成并下落，設置顆粒初始速度為2 m/s，0.25 s后全部花生莢果生成完畢，導管以0.05 m/s的速度向上運動，經(jīng)過5 s仿真時間后全部花生莢果顆粒在底部底板上靜止，形成花生莢果堆積角。

2.2 離散元仿真參數(shù)標定

2.2.1 確定顯著性影響參數(shù)

利用Design-Expert軟件進行Plackett-Burman試驗設計，試驗參數(shù)依照實際物理試驗測定結(jié)果，以實際花生莢果堆積角為目標值，通過Plackett-Burman試驗篩選出對目標值具有顯著性影響的參數(shù)[28]。分別將表2中的8個試驗參數(shù)的最大、最小值編碼為水平+1、-1。Plackett-Burman試驗方案及結(jié)果如表3所示。

表2 Plackett-Burman 試驗參數(shù)范圍表Tab. 2 Plackett-Burman test parameter range table

表3 Plackett-Burman試驗方案及結(jié)果Tab. 3 Plackett-Burman test protocol and results

在Design-Expert軟件中對試驗結(jié)果進行方差分析，得到仿真試驗參數(shù)顯著性結(jié)果，如表4所示。花生莢果—花生莢果滾動摩擦系數(shù)的P值<0.01，對仿真模擬試驗有極其顯著影響；花生莢果—花生莢果靜摩擦系數(shù)、花生莢果—鋼板靜摩擦系數(shù)的P值<0.05，對仿真模擬試驗有顯著影響；其他仿真試驗參數(shù)的P值>0.05，對仿真模擬試驗影響較小。

表4 Plackett-Burman試驗參數(shù)顯著性分析Tab. 4 Significance analysis of Plackett-Burmantest parameters

2.2.2 最陡爬坡試驗設計

在Plackett-Burman試驗的基礎上，對篩選出的3個顯著性因素(花生莢果—花生莢果靜摩擦系數(shù)、滾動摩擦系數(shù)，花生莢果—鋼板靜摩擦系數(shù))進行最陡爬坡試驗，以離散元仿真分析堆積角與實際堆積角的相對誤差作為評價標準，確定仿真參數(shù)最優(yōu)范圍。仿真試驗中，花生莢果泊松比取0.4，剪切模量取6.5 MPa；花生莢果—花生莢果碰撞恢復系數(shù)取0.275，花生莢果—鋼板碰撞恢復系數(shù)取0.275，花生莢果—鋼板滾動摩擦系數(shù)0.125。最陡爬坡試驗設計方案及結(jié)果如表5所示。仿真試驗結(jié)果顯示，3號試驗時，相對誤差最小，因此可以確定最優(yōu)參數(shù)范圍在3號附近，將3號作為中心點，2、4號分別當作低、高水平進行后續(xù)的Box-Behnken響應面試驗。

表5 最陡爬坡試驗設計方案及結(jié)果Tab. 5 Steepest climbing test design scheme and results

2.2.3 Box-Behnken試驗設計

以最陡爬坡試驗中3號當作中心點(0)，2號、4號分別當作低(-1)、高水平(+1)在Design-Expert軟件中進行顯著性參數(shù)的Box-Behnken試驗，設計方案及結(jié)果如表6所示。仿真試驗中其他非顯著性參數(shù)均按照最陡爬坡試驗中所用參數(shù)。

表6 Box-Behnken 試驗設計方案及結(jié)果Tab. 6 Box-Behnken experiment design scheme and results

在Design-Expert軟件中對Box-Behnken試驗結(jié)果進行二次響應面回歸分析，得到離散元仿真試驗堆積角二次回歸方程

θ=151.578-3 367.919Dx+214.719Ex+

4 273.733Gx+525DxEx-1 444.444DxGx-

344.444ExGx+2 523.056Dx2-771.389Ex2-

2 701.944 4Gx2

Box-Behnken試驗方差分析結(jié)果如表7所示，由分析結(jié)果可知Dx、Ex、DxGx、Dx2、Ex2、Gx2對堆積角影響極其顯著；Gx、DxEx對堆積角影響顯著；ExGx對堆積角影響不顯著。該回歸模型P<0.000 1，失擬項P=0.735 2>0.05，表明該回歸模型擬合度較好，無失擬現(xiàn)象發(fā)生?；貧w方程決定系數(shù)R2=0.987 5，校正決定系數(shù)R2=0.971 5，接近于1；回歸方程變異系數(shù)CV=1.18%。綜上分析數(shù)據(jù)所示，該回歸模型極其顯著，能夠反映真實情況，可用于目標堆積角的預測分析。

表7 Box-Behnken試驗回歸模型方差分析Tab. 7 Variation analysis of Box-Behnkentest quadratic model

2.2.4 仿真參數(shù)標定與試驗驗證

通過Design-Expert軟件中的優(yōu)化模塊，以實際試驗堆積角(31.63°)為目標值，對二次回歸方程式(8)進行優(yōu)化求解，從若干最優(yōu)參數(shù)組合中選擇一組與實際物理試驗測定數(shù)堆積相近的參數(shù)組合。花生莢果—花生莢果靜摩擦系數(shù)、花生莢果—花生莢果滾動摩擦系數(shù)、花生莢果—鋼板靜摩擦系數(shù)分別為0.74、0.24、0.58，其他非顯著性仿真參數(shù)取實際物理試驗所測定的平均值。

為驗證標定后花生莢果仿真參數(shù)的準確性，將上述標定后的物理參數(shù)作為離散元仿真參數(shù)，進行3次仿真模擬試驗驗證，得到花生莢果堆積角分別為32.29°、33.58°、31.75°，對仿真模擬試驗數(shù)據(jù)與實際試驗數(shù)據(jù)進行獨立樣本T檢驗，得到P=0.169>0.05，檢驗結(jié)果表明實際試驗堆積角與仿真參數(shù)標定后的仿真模擬試驗堆積角無顯著性差異，且實際試驗堆積角平均值31.63°與仿真模擬試驗堆積角平均值32.54°的相對誤差為2.877%，試驗對比如圖10所示，驗證了仿真模擬試驗的準確性。

(a) 實際物理試驗

同時對比了其他物料參數(shù)標定時所用試驗方法及標定結(jié)果，如表8所示。進一步表明標定后的參數(shù)可為相關研究提供參考。

表8 部分物料試驗方法及結(jié)果比較Tab. 8 Partial material test methods and comparison of results

3 結(jié)論

1) 通過實際試驗測得花生莢果基本物理參數(shù)，花生莢果外形尺寸為(27.575±2.049) mm×(14.39±0.955) mm×(13.646±0.914) mm，密度為465.08±26.51 kg/m，含水率為(19.28±0.335)%，容重為306.865±9.222 kg/m3，泊松比為0.4，彈性模量為17.46±5.72 MPa，剪切模量為6.5 MPa；采用HiSpec5型高速攝像機和自制摩擦系數(shù)測量儀測得花生莢果間的碰撞恢復系數(shù)、靜摩擦系數(shù)、滾動摩擦系數(shù)的平均值分別為0.287±0.043、0.823±0.041、0.283±0.043；花生莢果與Q235鋼板間的碰撞恢復系數(shù)、靜摩擦系數(shù)、滾動摩擦系數(shù)的平均值分別為0.275±0.064、0.583±0.015、0.104±0.024。

2) 以實際試驗測得的物理參數(shù)作為仿真模擬試驗參數(shù)選擇依據(jù)，開展Plackett-Burman試驗，篩選出對堆積角存在顯著性影響的因素為：花生莢果—花生莢果靜摩擦系數(shù)、滾動摩擦系數(shù)，花生莢果—鋼板靜摩擦系數(shù)，并進一步通過最陡爬坡試驗確定顯著性因素的取值范圍。

3) 開展Box-Behnken試驗，建立了堆積角與顯著性因素之間的二次回歸方程，并以實際物理試驗堆積角(31.63°)為目標值對方程進行求解，得到最佳仿真模擬參數(shù)：花生莢果—花生莢果靜摩擦系數(shù)、滾動摩擦系數(shù)、花生莢果—鋼板靜摩擦系數(shù)分別為0.74、0.24和0.58。

4) 對試驗分析后確定的最佳仿真參數(shù)進行仿真模擬試驗，對取得的仿真模擬值與實際試驗值進行獨立樣本T檢驗，得到P=0.169>0.05，表明實際試驗堆積角與仿真模擬試驗休止角無顯著性差異，且相對誤差為2.877%。驗證了仿真模擬試驗的準確性。通過對比其他物料參數(shù)標定時所用方法及試驗結(jié)果，進一步表明標定后的參數(shù)可為相關研究提供參考。