蔗段離散元仿真建模方法與參數(shù)標(biāo)定

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

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

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

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

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

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

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

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

隨機選取30根蔗段作為樣本，測量基本物理參數(shù)，取其平均值。采用卷尺、電子游標(biāo)卡尺(0.02 mm)、電子秤(0.01 g)測量蔗段幾何尺寸和質(zhì)量，采用排水法測量體積，采用電熱干燥法測量含水率。試驗測得蔗段幾何尺寸(長×直徑)的均值為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 蔗段物理堆積角測定

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

圖1 蔗段物理堆積試驗Fig. 1 Physical repose angle test of sugarcane segment

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 蔗段仿真堆積角測定

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

圖5 圖像處理Fig. 5 Image processing

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

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

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

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

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

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

仿真試驗時，需要用到蔗段和鋼板的基本物性參數(shù)和表面接觸參數(shù)，通過查閱文獻(xiàn)[33-39]對蔗段與鋼板力學(xué)特性的測定及離散元仿真中仿真參數(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ù)和滾動摩擦系數(shù)范圍分別為0.30～0.60、0.30～0.80 和 0.01～0.06；蔗段-鋼板的恢復(fù)系數(shù)、靜摩擦系數(shù)和滾動摩擦系數(shù)范圍分別為0.20～0.70、0.10～0.60 和 0～0.50。

3.2 Plackett-Burman(PB)試驗

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

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

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

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

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

3.3 最陡爬坡試驗

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

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

3.4 Box-Behnken(BB)試驗

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

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

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

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

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

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

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

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

圖7 堆積角仿真試驗驗證Fig. 7 Verification of repose angle by simulation test

4 結(jié)論

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

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

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