花生莢果離散元仿真參數標定*

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

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

0 引言

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

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

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

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

1 材料和方法

1.1 試驗材料與實驗儀器

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

1.2 物理參數測定

1.2.1 基本物理參數測定

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

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

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

1.2.2 泊松比

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

(1)

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

εX——寬度方向應變；

εY——長度方向應變；

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

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

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

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

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

(2)

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

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

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

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

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

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

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

(3)

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

1.3 接觸參數測定

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

1.3.1 碰撞恢復系數測定

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

(4)

式中：e——物體的恢復系數；

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

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

(5)

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

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

圖2 恢復系數測試原理圖

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

圖3 跌落碰撞試驗裝置

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

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

1.3.2 靜摩擦系數測定

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

圖4 摩擦系數測量原理示意

(6)

μ1=f/N=tanθ

(7)

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

μ1——花生莢果靜摩擦系數；

m——花生莢果質量，kg；

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

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

圖5 摩擦系數斜面儀

1.3.3 動摩擦系數測定

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

1.3.4 堆積角測定

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

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

圖7 單側堆積角邊界擬合

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

2.1 離散元仿真模型建立

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

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

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

圖9 堆積角測定裝置

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

2.2 離散元仿真參數標定

2.2.1 確定顯著性影響參數

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

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

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

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

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

2.2.2 最陡爬坡試驗設計

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

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

2.2.3 Box-Behnken試驗設計

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

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

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

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

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

2.2.4 仿真參數標定與試驗驗證

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

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

(a) 實際物理試驗

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

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

3 結論

1) 通過實際試驗測得花生莢果基本物理參數，花生莢果外形尺寸為(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型高速攝像機和自制摩擦系數測量儀測得花生莢果間的碰撞恢復系數、靜摩擦系數、滾動摩擦系數的平均值分別為0.287±0.043、0.823±0.041、0.283±0.043；花生莢果與Q235鋼板間的碰撞恢復系數、靜摩擦系數、滾動摩擦系數的平均值分別為0.275±0.064、0.583±0.015、0.104±0.024。

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

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

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