任甲輝，武 濤，劉慶庭，莫王玉潔，李 珂

(華南農(nóng)業(yè)大學(xué) 工程學(xué)院,廣東 廣州 510642)

甘蔗是我國(guó)重要的糖料作物，2019年，我國(guó)糖料蔗種植面積約153.33萬(wàn)hm2，位居世界第三[1]。目前甘蔗機(jī)械化收獲技術(shù)存在含雜率高、損失大的突出問(wèn)題，難以得到糖廠和農(nóng)戶(hù)的認(rèn)可，嚴(yán)重制約了我國(guó)甘蔗收獲機(jī)的推廣[2-3]。甘蔗在收獲過(guò)程中，物料間、物料與設(shè)備間的相互作用十分復(fù)雜，采用離散元法研究顆粒間的相互作用及運(yùn)動(dòng)規(guī)律，有利于揭示顆粒間相互作用機(jī)理，進(jìn)而改善甘蔗輸送、除雜等關(guān)鍵環(huán)節(jié)的作業(yè)效率，對(duì)提高甘蔗收獲機(jī)性能具有重要意義。顆粒仿真參數(shù)的標(biāo)定，有利于提高離散元仿真模型的精度，從而提高離散元法在顆粒間相互作用機(jī)理及運(yùn)動(dòng)規(guī)律研究中的準(zhǔn)確度[4-5]。

離散元法是研究物料運(yùn)動(dòng)的有效工具，近年來(lái)，在農(nóng)業(yè)工程研究中呈現(xiàn)出良好的應(yīng)用前景[6]。劉凡一等[7]以小麥顆粒物理休止角為響應(yīng)值，基于響應(yīng)面優(yōu)化標(biāo)定了小麥離散元仿真參數(shù)。侯占豐等[8]通過(guò)物理試驗(yàn)的方法測(cè)定了冰草種子的基本物性參數(shù)和接觸參數(shù)，以物理休止角為最優(yōu)目標(biāo)值，對(duì)離散元仿真參數(shù)進(jìn)行尋優(yōu)，獲得最佳參數(shù)組合。吳孟宸等[9]對(duì)花生種子顆粒堆積過(guò)程進(jìn)行仿真，利用堆積角試驗(yàn)對(duì)種間摩擦因數(shù)及種子與材料間滾動(dòng)摩擦因數(shù)進(jìn)行離散元仿真標(biāo)定，通過(guò)臺(tái)架試驗(yàn)驗(yàn)證花生模型及其仿真參數(shù)的可靠性。武濤等[10]采用顆粒間黏結(jié)力的“Hertz-Mindlin with JKR”接觸模型，進(jìn)行土壤堆積角仿真試驗(yàn)，確定了樣品土壤的離散元接觸參數(shù)與接觸模型參數(shù)。王黎明等[11]測(cè)定了不同含水率下豬糞的堆積角，通過(guò)物理堆積試驗(yàn)與仿真方法對(duì)豬糞接觸參數(shù)進(jìn)行了標(biāo)定。劉文政等[12]利用物理試驗(yàn)與仿真試驗(yàn)相結(jié)合的方法對(duì)微型馬鈴薯的離散元仿真參數(shù)進(jìn)行了優(yōu)化和標(biāo)定。國(guó)內(nèi)外學(xué)者還對(duì)水稻[13-14]、大豆[15]、玉米[16-17]等常見(jiàn)谷物以及葡萄[18]、蘋(píng)果[19-20]等水果進(jìn)行了參數(shù)標(biāo)定。綜上所述，目前農(nóng)業(yè)物料離散元仿真模型參數(shù)標(biāo)定主要集中在糧油種子、土壤、肥料以及生物質(zhì)顆粒等球形/類(lèi)球形微小顆粒，以及水果等類(lèi)球形大顆粒。堆積角試驗(yàn)是進(jìn)行離散元仿真參數(shù)標(biāo)定的主要方法，對(duì)于上述顆粒，對(duì)堆積角影響顯著性較大的參數(shù)主要為恢復(fù)系數(shù)、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)。切段式甘蔗收獲機(jī)的收獲物為甘蔗段，與上述傳統(tǒng)球形/類(lèi)球形物料不同，蔗段屬于大長(zhǎng)徑比稈狀物料，對(duì)這類(lèi)物料堆積角影響較顯著的仿真參數(shù)尚不清楚，影響了離散元法在蔗段運(yùn)動(dòng)學(xué)分析和CFDEDM耦合除雜機(jī)理研究方面的應(yīng)用。因此，亟需進(jìn)行蔗段離散元建模方法和參數(shù)標(biāo)定研究。

本文以蔗段為研究對(duì)象，選用多球聚合模型和XML的方法構(gòu)建蔗段仿真模型，根據(jù)國(guó)內(nèi)外研究確定蔗段仿真參數(shù)范圍，以物理試驗(yàn)獲取的蔗段堆積角為響應(yīng)值，通過(guò)Plackett-Burman試驗(yàn)、最陡爬坡試驗(yàn)和Box-Behnken試驗(yàn)對(duì)離散元仿真參數(shù)進(jìn)行標(biāo)定，得到最優(yōu)參數(shù)組合，最終通過(guò)物理試驗(yàn)與仿真試驗(yàn)的對(duì)比驗(yàn)證，確定離散元仿真參數(shù)的可靠性。以期為蔗段等大長(zhǎng)徑比非球形物料離散元仿真參數(shù)標(biāo)定提供參考。

1 蔗段基本物性參數(shù)與堆積角測(cè)定

1.1 蔗段基本物性參數(shù)測(cè)定

選取廣東雷州半島南光6隊(duì)農(nóng)場(chǎng)種植的‘桂糖49’甘蔗為研究對(duì)象，采用華南農(nóng)業(yè)大學(xué)自主研制的HN4GDL-194型切段式甘蔗聯(lián)合收獲機(jī)為收獲設(shè)備，對(duì)切段后的甘蔗進(jìn)行基本物性參數(shù)(幾何尺寸、密度、含水率)測(cè)量。

隨機(jī)選取30根蔗段作為樣本，測(cè)量基本物理參數(shù)，取其平均值。采用卷尺、電子游標(biāo)卡尺(0.02 mm)、電子秤(0.01 g)測(cè)量蔗段幾何尺寸和質(zhì)量，采用排水法測(cè)量體積，采用電熱干燥法測(cè)量含水率。試驗(yàn)測(cè)得蔗段幾何尺寸(長(zhǎng)×直徑)的均值為283.08 mm×32.06 mm，標(biāo)準(zhǔn)差分別為2.23和1.60，蔗段的平均質(zhì)量為 192 g，密度為 1.19 g/cm3，含水率（w）為 73.15%。

1.2 蔗段物理堆積角測(cè)定

堆積角是表征散裝物料流動(dòng)、摩擦等性質(zhì)的宏觀參數(shù)，與物料種類(lèi)、長(zhǎng)徑比、形狀和含水率等因素有關(guān)，因此堆積角試驗(yàn)常被用作散裝物料的離散元標(biāo)定[10]。蔗段物理堆積角的測(cè)量方法為：采用華南農(nóng)業(yè)大學(xué)自主研制的甘蔗田間運(yùn)輸車(chē)，將蔗段從田間運(yùn)輸車(chē)的輸送臂卸到公路運(yùn)輸車(chē)內(nèi)，蔗段從相對(duì)高度1 500 mm處落下，待所有蔗段靜止后，在公路運(yùn)輸車(chē)內(nèi)形成穩(wěn)定的甘蔗堆，用角度測(cè)量?jī)x測(cè)量蔗段堆積角，試驗(yàn)重復(fù)10次，獲得蔗段物理試驗(yàn)堆積角為39.40°～45.80°，平均值為42.70°，標(biāo)準(zhǔn)差為2.6。蔗段物理堆積試驗(yàn)，如圖1所示。

圖1 蔗段物理堆積試驗(yàn)Fig. 1 Physical repose angle test of sugarcane segment

2 蔗段堆積角離散元仿真模型創(chuàng)建

2.1 非球形顆粒模型構(gòu)建方法

球形顆粒形狀簡(jiǎn)單、規(guī)則，只有1個(gè)尺寸參數(shù)，因此在離散元仿真軟件中，一般選用球形顆粒建立物料的仿真模型?？紤]到離散元仿真效率，通常球形/類(lèi)球形物料離散元模型可以用單個(gè)球形顆粒進(jìn)行簡(jiǎn)化[21]。但實(shí)際中，大多數(shù)物料都是不規(guī)則的，無(wú)法通過(guò)單一球形進(jìn)行簡(jiǎn)化建模，為使離散元仿真模型更符合實(shí)際物料特性，需采用非球形顆粒模型的構(gòu)建方法。

非球形顆粒建模通常采用多球聚合模型(Multisphere method，MSM)[22-24]和多球黏結(jié)模型 (Bonded particle method，BPM)[25-26]。多球聚合模型 (圖2a)由多個(gè)顆粒堆疊聚合而成，該離散元模型與實(shí)際物料輪廓相對(duì)吻合，且在EDEM軟件中被視為一個(gè)獨(dú)立體。多球黏結(jié)模型(圖2b)由多個(gè)直徑相同的顆粒通過(guò)“黏結(jié)鍵”粘連而成，顆粒間相互獨(dú)立，半徑越小，越接近實(shí)際物料輪廓，但是會(huì)增加仿真時(shí)間，適用于DEM-CFD耦合仿真[13,27]。

圖2 非球形顆粒建模方法Fig. 2 Modeling method of non-spherical particles

非球形顆粒模型填充方法有手動(dòng)填充、API顆粒替換和XML等方法。手動(dòng)填充常用于少量球形顆粒填充的簡(jiǎn)單模型。API顆粒替換和XML方法常用于大量球形顆粒填充建模，可提高構(gòu)建顆粒模型的速度和精度。API顆粒替換方法需要調(diào)用相應(yīng)插件完成小顆粒替換大顆粒并建立黏結(jié)鍵的過(guò)程，小顆粒之間相互獨(dú)立，常用于磨損、破裂方面的研究[28]。XML方法是在EDEM后處理模塊導(dǎo)出.xml 格式文件中，調(diào)整數(shù)據(jù)文件中顆粒的坐標(biāo)得到目標(biāo)顆粒，所有顆粒被視為一個(gè)整體[29]，常用于不考慮磨損、破裂等問(wèn)題的模型構(gòu)建。

2.2 蔗段模型構(gòu)建

蔗段屬于稈狀物料模型，在簡(jiǎn)化此類(lèi)物料模型結(jié)構(gòu)時(shí)，不考慮物料的外觀細(xì)微變化[30]。根據(jù)前期物理試驗(yàn)得到的蔗段基本物理參數(shù)，采用Solidworks建立幾何模型，將蔗段模型保存為.stl格式導(dǎo)入EDEM中。本文選用多球聚合模型和XML的方法構(gòu)建蔗段仿真模型，定義顆粒工廠在蔗段內(nèi)部上端，動(dòng)態(tài)生成的顆粒自由下落至模型內(nèi)部，顆粒間采用Hertz-Mindin無(wú)滑移接觸模型。利用154個(gè)半徑為8 mm的顆粒自動(dòng)填充蔗段模型，填充完成后，導(dǎo)出所有球形顆粒的中心坐標(biāo).csv格式文件，對(duì)文件中的數(shù)據(jù)進(jìn)行處理。將仿真時(shí)間歸零，導(dǎo)出格式為.xml的DEC文件，將EDEM中新建的paricle項(xiàng)目導(dǎo)入.xml格式文件后，軟件自動(dòng)生成的蔗段離散元仿真模型，可作為顆粒模板，進(jìn)行后續(xù)仿真調(diào)用，蔗段的離散元模型如圖3所示。

圖3 蔗段的仿真模型Fig. 3 Simulation model of sugarcane segment

2.3 蔗段堆積角仿真模型構(gòu)建

根據(jù)蔗段的物理堆積角試驗(yàn)，簡(jiǎn)化仿真試驗(yàn)幾何模型。物理試驗(yàn)時(shí)，輸送臂只起到將蔗段提升至一定高度的作用，蔗段從輸送帶上輸送、拋散、下落，散開(kāi)后豎直下落的最大圓周直徑約400 mm，豎直下落的初始位置與輸送臂末端高度差約300 mm，與運(yùn)輸車(chē)箱底板高度差約1 500 mm，蔗段下落到運(yùn)輸車(chē)內(nèi)，穩(wěn)定后蔗段堆的最大散落直徑約1 500 mm。仿真時(shí)，為了提高效率，簡(jiǎn)化去掉輸送臂和運(yùn)輸車(chē)，將蔗段豎直下落的初始位置設(shè)置為漏斗下開(kāi)口，豎直下落的初始位置正上方300 mm處設(shè)置為顆粒工廠(漏斗上開(kāi)口)，公路運(yùn)輸車(chē)簡(jiǎn)化為單一接料板。

在EDEM中建立漏斗(上、下開(kāi)口直徑分別為500 mm 和 400 mm、高度 300 mm、下端面離地高度為 1 500 mm)和接料板 (直徑 1 500 mm)，設(shè)置顆粒工廠在漏斗上方口徑處，總生成量為1 400根，生成速率為100根/s，仿真時(shí)間為16 s，仿真蔗段的堆積過(guò)程，如圖4所示。

圖4 堆積角仿真模型Fig. 4 Repose angle of simulation model

2.4 蔗段仿真堆積角測(cè)定

蔗段仿真堆積角圖像處理過(guò)程，如圖5所示。

圖5 圖像處理Fig. 5 Image processing

仿真試驗(yàn)中，為了更準(zhǔn)確地測(cè)量仿真堆積角，采用計(jì)算機(jī)圖像處理技術(shù)對(duì)其進(jìn)行處理。首先在EDEM中截取仿真蔗段堆的正視圖，如圖5a所示，再利用Matlab軟件對(duì)獲得的圖像進(jìn)行灰度化、二值化（圖5b）、邊界像素點(diǎn)的搜索與提取，最后通過(guò)線性擬合獲取擬合方程，得到方程斜率K[31]。蔗段堆積角 (θ)計(jì)算公式[32]為：

式中：K為擬合曲線方程的斜率。

由于蔗段堆X、Y方向兩側(cè)傾角存在一定差異，因此測(cè)量視圖兩側(cè)的傾角，取其平均值為仿真堆積角，如圖6所示。

圖6 蔗段離散元仿真模型堆Fig. 6 DEM simulation model stack of sugarcane segment

3 蔗段離散元仿真參數(shù)標(biāo)定

3.1 仿真參數(shù)范圍確定

仿真試驗(yàn)時(shí)，需要用到蔗段和鋼板的基本物性參數(shù)和表面接觸參數(shù)，通過(guò)查閱文獻(xiàn)[33-39]對(duì)蔗段與鋼板力學(xué)特性的測(cè)定及離散元仿真中仿真參數(shù)的設(shè)置，確定本研究中各仿真參數(shù)的取值范圍：蔗段的泊松比為 0.30～0.38、剪切模量為 3.0～18.6 MPa、密度為1.19 g/cm3；鋼板的泊松比為0.30、剪切模量為 7.9×104MPa、密度為 7.85 g/cm3；蔗段-蔗段的恢復(fù)系數(shù)、靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)范圍分別為0.30～0.60、0.30～0.80 和 0.01～0.06；蔗段-鋼板的恢復(fù)系數(shù)、靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)范圍分別為0.20～0.70、0.10～0.60 和 0～0.50。

3.2 Plackett-Burman(PB)試驗(yàn)

為了篩選出對(duì)蔗段堆積角存在顯著性影響的仿真參數(shù)，開(kāi)展PB試驗(yàn)。利用Design-Expert.11軟件對(duì)PB試驗(yàn)進(jìn)行設(shè)計(jì)，以蔗段物理試驗(yàn)堆積角為響應(yīng)值，進(jìn)行顯著性參數(shù)篩選。選取表1中的8個(gè)仿真參數(shù)的最大值、最小值作為高、低2個(gè)水平，分別編碼為水平1、-1，如表1所示。以中間水平0為中心點(diǎn)，進(jìn)行13組試驗(yàn)。

表1 Plackett-Burman試驗(yàn)參數(shù)列表Table 1 List of Plackett-Burman test parameters

PB試驗(yàn)方案及結(jié)果如表2所示，利用Design-Expert.11軟件對(duì)試驗(yàn)結(jié)果進(jìn)行方差分析，得到各仿真參數(shù)顯著性結(jié)果如表3所示。由表3可知，蔗段-蔗段靜摩擦系數(shù)(X4)和蔗段-蔗段滾動(dòng)摩擦系數(shù)(X5)對(duì)仿真堆積角的影響極顯著(P＜0.01)；蔗段泊松比(X1)對(duì)仿真堆積角的影響顯著(P＜0.05)；蔗段剪切模量、蔗段-蔗段恢復(fù)系數(shù)、蔗段-鋼恢復(fù)系數(shù)、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)對(duì)仿真堆積角的影響不顯著。因此后續(xù)試驗(yàn)設(shè)計(jì)中只考慮3個(gè)影響顯著的因素，對(duì)仿真堆積角影響不顯著的參數(shù)取PB試驗(yàn)的中間水平。

表2 Plackett-Burman試驗(yàn)方案及結(jié)果1)Table 2 Plackett-Burman test scheme and results

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

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

為快速確定各仿真參數(shù)最優(yōu)值所在區(qū)域范圍，對(duì)PB試驗(yàn)篩選出的3個(gè)顯著性參數(shù)進(jìn)行爬坡試驗(yàn)，包括蔗段的泊松比、蔗段-蔗段的靜摩擦系數(shù)、蔗段-蔗段的滾動(dòng)摩擦系數(shù)。試驗(yàn)以仿真堆積角與物理堆積角的相對(duì)誤差作為評(píng)價(jià)指標(biāo)，相對(duì)誤差最小的組可確定為各仿真參數(shù)最優(yōu)值區(qū)域間的中心點(diǎn)。最陡爬坡試驗(yàn)方案設(shè)計(jì)及結(jié)果（表4）表明，隨著3個(gè)顯著性參數(shù)數(shù)值的增加，仿真堆積角與物理試驗(yàn)堆積角的相對(duì)誤差先減小后增大，在4號(hào)水平時(shí)，相對(duì)誤差最小為0.62%，可以確定最優(yōu)值區(qū)域在4號(hào)水平附近，因此取4號(hào)作為中心點(diǎn)，3號(hào)、5號(hào)分別作為低、高水平進(jìn)行后續(xù)響應(yīng)面設(shè)計(jì)。

表4 最陡爬坡試驗(yàn)方案設(shè)計(jì)及結(jié)果Table 4 Design scheme and results of the steepest ascent search test

3.4 Box-Behnken(BB)試驗(yàn)

根據(jù)最陡爬坡試驗(yàn)結(jié)果，以表4中4號(hào)作為中心點(diǎn)(0)，3號(hào)和5號(hào)分別為低水平(-1)和高水平(1)進(jìn)行顯著性參數(shù)Box-Behnken試驗(yàn)，設(shè)計(jì)方案及結(jié)果如表5所示。

表5 Box-Behnken試驗(yàn)設(shè)計(jì)方案及結(jié)果1)Table 5 Box-Behnken test design scheme and results

利用Design-Expert.11件對(duì)BB試驗(yàn)結(jié)果進(jìn)行二次多項(xiàng)式回歸擬合，得到蔗段仿真堆積角(θ)與3個(gè)顯著性參數(shù)的回歸方程：

該回歸模型方差分析結(jié)果如表6所示，由表6分析可知，蔗段泊松比(X1)、蔗段-蔗段靜摩擦系數(shù)(X4)、蔗段-蔗段滾動(dòng)摩擦系數(shù)二次項(xiàng)對(duì)試驗(yàn)堆積角結(jié)果影響極其顯著，蔗段-蔗段靜摩擦系數(shù)二次項(xiàng)對(duì)試驗(yàn)堆積角結(jié)果影響顯著，其他各項(xiàng)影響均不顯著。該堆積角擬合模型P＜0.01，說(shuō)明該模型擬合度極顯著；失擬項(xiàng)P=0.318 3＞0.05，說(shuō)明模型擬合良好，無(wú)失擬因素存在。該回歸模型的決定系數(shù)為 0.968 0，校正決定系數(shù)為 0.910 4，兩者均接近1，說(shuō)明回歸方程的預(yù)測(cè)值與實(shí)際值擬合良好；變異系數(shù)為1.41%，說(shuō)明試驗(yàn)較可靠；精確度為15.186 6，說(shuō)明模型具有良好精度。綜上，該回歸模型能真實(shí)、可靠地反應(yīng)實(shí)際情況，可用來(lái)進(jìn)一步預(yù)測(cè)顆粒堆的堆積角。

表6 Box-Behnken試驗(yàn)回歸模型方差分析Table 6 Analysis of variance of the Box-Behnken test regression model

3.5 最優(yōu)參數(shù)組合與試驗(yàn)驗(yàn)證

利用Design-Expert.11軟件的優(yōu)化模塊，以物理堆積角42.70°為目標(biāo)值，對(duì)回歸方程進(jìn)行優(yōu)化求解，得到多組最佳參數(shù)。從多組最佳參數(shù)組合中選取1組與物理堆積角平均值最接近的參數(shù)組合。最終選取的最佳優(yōu)化組合為：蔗段泊松比0.35、蔗段-蔗段靜摩擦系數(shù)0.53、蔗段-蔗段滾動(dòng)摩擦系數(shù)0.04，其余非顯著性參數(shù)取中間水平。

為了驗(yàn)證蔗段離散元標(biāo)定參數(shù)的準(zhǔn)確性，采用上述最優(yōu)參數(shù)組合進(jìn)行蔗段堆積角仿真試驗(yàn)，重復(fù)3次仿真模擬，得到的蔗段仿真堆積角分別為42.69°、43.06°和43.65°。對(duì)仿真結(jié)果與物理試驗(yàn)值進(jìn)行獨(dú)立樣本t檢測(cè)，得到P=0.196＞0.05，表明標(biāo)定后的仿真堆積角與物理試驗(yàn)堆積角無(wú)顯著性差異，仿真試驗(yàn)均值(43.13°)與物理試驗(yàn)均值(42.70°)相對(duì)誤差為0.99%，進(jìn)一步驗(yàn)證了蔗段離散元標(biāo)定參數(shù)的可靠性。蔗段堆積角仿真試驗(yàn)驗(yàn)證如圖7所示。

圖7 堆積角仿真試驗(yàn)驗(yàn)證Fig. 7 Verification of repose angle by simulation test

4 結(jié)論

研究通過(guò)物理試驗(yàn)和參考國(guó)內(nèi)外文獻(xiàn)，確定了蔗段基本參數(shù)及取值范圍；并通過(guò)Plackett-Burman試驗(yàn)，對(duì)仿真參數(shù)進(jìn)行顯著性篩選。結(jié)果表明：蔗段泊松比、蔗段-蔗段靜摩擦系數(shù)和蔗段-蔗段滾動(dòng)摩擦系數(shù)對(duì)蔗段堆積角影響顯著。與以往研究[13-20]相比較，發(fā)現(xiàn)泊松比對(duì)大長(zhǎng)徑比稈狀物料堆積角的影響比球類(lèi)、類(lèi)球形物料更為顯著。

通過(guò)Box-Behnken試驗(yàn)，以物理試驗(yàn)堆積角42.70°為最佳目標(biāo)值，得到最優(yōu)參數(shù)組合為：蔗段泊松比0.35、蔗段-蔗段靜摩擦系數(shù)0.53、蔗段-蔗段滾動(dòng)摩擦系數(shù)0.04。

通過(guò)獨(dú)立樣本t檢驗(yàn)，對(duì)最優(yōu)參數(shù)組合進(jìn)行驗(yàn)證。結(jié)果表明仿真堆積角與物理試驗(yàn)堆積角無(wú)顯著性差異，兩者相對(duì)誤差為0.99%，進(jìn)一步驗(yàn)證了蔗段離散元標(biāo)定參數(shù)的可靠性。