聯合收獲油菜脫出物離散元仿真參數標定與試驗

舒彩霞，楊 佳，萬星宇，袁佳誠，廖宜濤，廖慶喜

（1.華中農業大學工學院，武漢 430070；2.農業農村部長江中下游農業裝備重點實驗室，武漢 430070）

0 引 言

油菜是中國主要油料作物，其生產種植面積占全國油菜總種植面積85%以上，全程全面高質高效推進油菜機械化發展是保障食用油安全的重要措施。油菜聯合收獲是油菜機械化收獲的主要方式之一，可一次性完成切割、脫粒、分離、清選等工序，直接獲得干凈的油菜籽粒，生產效率高。在收獲時，過早收獲會產生脫粒不凈、青籽多等問題；過晚收獲容易造成落粒、損失率增加；因此，最適宜收獲時期是田間油菜呈黃熟期后至完熟期之間，90%以上的油菜角果變成黃色和褐色，主分支向上收攏，果莢中黑籽粒在75%左右。

為提高油菜聯合收獲機清選裝置的清選性能，國內外學者對清選裝置進行了研究。油菜籽粒在進行清選作業時，脫出物中的籽粒、短莖稈及莢殼間的相互作用關系難以通過數學模型進行分析，采用離散元法（Discrete Element Method，DEM）和計算機流體力學（Computational Fluid Dynamics，CFD）對關鍵部件進行數值模擬可以有效模擬出物料運動軌跡，替代繁雜的臺架試驗，省時、省力，縮短研發周期，降低成本。在進行油菜脫出物清選仿真試驗時，脫出物間精確的離散元接觸參數和特征參數可以提高氣固耦合清選仿真中顆粒間相互作用及運動規律的準確度，進而實現關鍵部件的參數優化。

為提高仿真模擬的準確性，使仿真結果更加接近實際作業過程，國內外學者基于堆積角試驗針對不同物料的離散元接觸參數進行了標定。Guo等通過物理試驗與仿真試驗相結合方法建立了香蕉莖稈力學模型并對其生物力學特性進行了研究。馬彥華等采用6個接觸參數測定試驗，以此確定優化試驗的參數區間，依次通過Plackett-Burman試驗、最陡爬坡試驗以及Box-Behnken 試驗對苜蓿秸稈的離散元仿真參數進行標定。侯占峰等通過物理試驗和仿真試驗相結合的方法對冰草種子進行了參數標定。石林榕等通過堆積試驗、查閱參考文獻和計算了馬鈴薯基本物理參數和接觸力學參數。綜上所述，國內外學者為不同作物的數值模擬提供了可參考的離散元參數。由于油菜脫出物組分糅雜、成熟度不一致、物料間特性差異大，國內外學者圍繞油菜脫出物的空氣動力學特性、特征參數及組分配比進行了研究。陳立等對油菜脫出物中籽粒、莖稈、莢殼及輕雜的懸浮速度進行了測定。Zhan等對整株植物的抗裂莢性、抗主莖剪斷性和共振頻率等生物力學特性進行研究。雷小龍針對油菜籽粒的特征參數及接觸參數進行了大量研究。袁佳誠等開展了脫出物質量比試驗，發現隨著切碎滾筒轉速增加，脫出物中油菜莖稈平均長度逐漸降低，增加了脫出物中短莖稈的質量比，影響旋風分離清選裝置的清選性能。相關研究可為仿真模擬的基礎參數提供依據，但相關研究主要是對油菜單一組分的部分物理特性和油菜籽粒的部分特征參數及接觸參數進行了研究，對于油菜脫出物（籽粒、短莖稈、莢殼）離散元接觸模型參數鮮有報道。

本文采用物理試驗與仿真試驗相結合的方法進行相關參數標定，對油菜脫出物組分間的靜摩擦系數，碰撞恢復系數進行標定，采用堆積試驗對莖稈特征參數、莖稈間接觸參數及莖稈-鋼板接觸參數進行標定，在優化標定后參數下，通過DEM-CFD氣固耦合試驗與臺架試驗驗證了標定參數準確性，以期為油菜聯合收獲機械化仿真模擬提供離散元模型基礎參數。

1 聯合收獲油菜脫出物物料特性測量與分析

1.1 試驗材料

2021年5月于華中農業大學現代農業科技試驗基地中采用自主研發的油菜聯合收獲機開展田間試驗，試驗材料為機直播華油雜62號油菜，平均種植密度為35～40 株/m，平均株高為1 655.44 mm，平均分支數為6。為準確反應莖稈的實際機械物理特性，采用五點取樣法獲得油菜莖稈，分別在1 m的測試區隨機人工選取油菜植株，從離地350 mm開始，每隔100 mm截取油菜莖稈；收集脫出物（如圖1所示）中莢殼及籽粒，將收集的材料分選后作為試驗材料，按照其種類人工分選后并稱量，其中籽粒、短莖稈和莢殼的質量比為2∶3∶5，含水率分別為30.41%～37.23%、40.92%～48.32%和66.67%～73.53%。

圖1 油菜脫出物 Fig.1 The rape threshing mixture

1.2 莖稈特征參數測量與分析

莖稈材料特征參數包括密度、泊松比、彈性模量和剪切模量等。通過MNT951221型數顯游標卡尺（測量誤差為0.02 mm）對收集的油菜莖稈直徑進行測量，測得油菜莖稈的平均直徑為10.05 mm。根據文獻[23]，通過體積計算和質量測定計算得到油菜莖稈密度為494 kg/m。

通過質構儀（TMS-Pro型，美國FTC公司）對隨機選取的油菜莖稈進行徑向單軸平板全壓縮試驗，試驗采用直徑75 mm 的平板壓頭，加載速度為50 mm/s，如圖2a所示。莖稈出現破裂聲時停止加載，采用定義法測量油菜莖稈泊松比；利用數顯游標卡尺記錄被擠壓方向及與其垂直方向的直徑變形量，使用質構儀記錄載荷-位移數據（如圖2b所示）。重復10次，通過式（1）～（3）計算得到油菜莖稈的彈性模量、泊松比和剪切模量平均值分別為4.44 MPa、0.45、1.53 MPa。

圖2 油菜莖稈特征參數測量 Fig.2 Measurement of feature parameters of rape stalk

式中為物料的彈性模量，MPa；為最大壓應力，Pa；為線應變；F為最大載荷，N；、為莖稈外徑和內徑，mm；為試樣接觸長度，75mm；Δ為莖稈直徑變化量，mm；Δ為莖稈壓縮前直徑，mm。

式中為物料的泊松比；為被擠壓方向應變；為與被擠壓垂直方向應變；為油菜莖稈試驗前被擠壓方向初始直徑，mm；Δ為油菜莖稈試驗后被擠壓方向直徑變化量，mm；為油菜莖稈試驗前與被擠壓垂直方向的初始直徑，mm；Δ為油菜莖稈試驗后與被擠壓垂直方向的直徑變化量，mm。

式中為物料的剪切模量，MPa。

1.3 脫出物接觸參數測量與分析

脫出物接觸參數主要包括靜摩擦系數、動摩擦系數及碰撞恢復系數等。油菜籽粒細小且脫出物生物學特性差異明顯、形態各異。參數參考文獻[14]，本文通過斜面法、自由跌落試驗對脫出物組分間靜摩擦系數和碰撞恢復系數進行測量，采用堆積角逼近的方式對莖稈與莖稈、鋼板動摩擦系數進行標定，明確物料組分間的接觸參數范圍。

采用斜面法測量不同物料間的靜摩擦系數。采用摩擦系數測定裝置，待測板為300 mm×300 mm× 3 mm鋼板，試驗裝置如圖3所示。當物體靜止于斜面上時：

圖3 物料間靜摩擦系數測量試驗 Fig.3 Static friction coefficient test between materials

式中為物料的重力，N；為斜面與水平面的夾角，(°)；為待測物料與物料種群板的靜摩擦系數。

通過不斷增加斜面傾角，當物體剛好有滑動趨勢時有：

此時斜面傾角正切值便是此時接觸物體間的靜摩擦系數。

在測定物料間靜摩擦系數時，用鑷子將油菜籽粒、莢殼、短莖稈等分別放至于物料種群板上，緩慢提高物料種群板，當物料出現向下運動時停止抬升物料種群板，記錄此時的角度，通過式（5），即可得到物料間的靜摩擦系數，每組試驗重復20次。根據文獻[14]，籽粒間靜摩擦系數取0.5。在測量時，為了充分模擬物料與接觸材料的接觸特性，使試驗數據有良好的適應性，物料每次放置位置盡量不重復，得到的油菜脫出物組分間的靜摩擦系數如表1所示。

表1 接觸物料間的靜摩擦系數 Table 1 Static friction coefficient between contact materials

采用自由跌落試驗測量不同油菜脫出物間碰撞恢復系數。將所測物料均簡化為質點，物料從高度開始掉落，與地面接觸板接觸時（碰撞開始）質心速度為，受到固定接觸面的碰撞沖量作用，物料變形增大，速度逐漸減小直至為0；此后彈性變形開始恢復，物料開始反向上升，物料離開接觸面時的瞬間獲得反向速度，并一直上升至高度處速度為0（碰撞結束）。碰撞過程如圖4所示，由沖量定理可以求出碰撞恢復系數表達式，如式（6）所示。

圖4 碰撞恢復系數測量原理圖 Fig.4 Schematic diagram of collision recovery coefficient test

式中為物料質量，g；、為彈性變形前后的碰撞沖量，N·s；為碰撞恢復系數。

為清楚觀察到物料的彈起高度，經過多次重復試驗，隨機選取莖稈制成20 mm的圓柱顆粒進行自由跌落試驗，物料自地面接觸板正上方200 mm處自由落體；分別以莖稈種群板、莢殼種群板以及鋼板作為地面接觸板進行測量。將網格紙粘于鋼板上作為背景記錄物料下降高度及彈起高度，采用美國VRI的Phantom系列高速攝影機正視放置于網格紙正前方，用鑷子夾持物料自200 mm處釋放，每組試驗重復20次，將拍攝的照片導入電腦，通過式（6）計算出不同物料的碰撞恢復系數，如表2所示，根據文獻[14]籽粒間碰撞恢復系數取0.6。

表2 接觸物料間的碰撞恢復系數 Table 2 Collision recovery coefficient between contact materials

2 油菜聯合收獲莖稈顆粒離散元參數標定

由于油菜莖稈為各向異性材料，因此本文采用仿真逼近的方法對油菜莖稈接觸參數進行標定。堆積角反應了顆粒物料流動及摩擦特性，與接觸材料及自身物理特性相關。因此，堆積角物理試驗常常被用來顆粒離散元參數標定。利用物理試驗標定得到的參數范圍開展油菜莖稈堆積試驗，不斷調整油菜莖稈間接觸參數使仿真堆積角逼近實際堆積角，對應得到油菜莖稈間準確接觸參數。

2.1 油菜莖稈實際堆積角測量

采用圓筒提升法獲得油菜莖稈堆來測量堆積角。參考文獻[32]，將收集的莖稈制成長度為20、30、40，50和60 mm莖稈顆料，經過多次重復試驗，每種長度的莖稈各150 粒，共750 粒。試驗所用的圓筒為內徑200 mm、高度300 mm的鋼制圓筒，將試樣填滿圓筒并放置于萬能材料試驗機平臺，將圓筒以0.05 m/s勻速提升，試驗重復3次。

正視拍照后將物理試驗圖片進行處理，應用Matlab中數字圖像處理軟件讀取顆粒堆邊緣圖像，對圖像進行灰度化、二值化，輸出二值圖，提取二值圖邊界輪廓，掃描輪廓圖像上每一個像素，記錄白點的坐標及個數，通過最小二乘法對所記錄的白點進行線性擬合，在Matlab讀取擬合直線的斜率，求得顆粒堆邊緣輪廓的傾角。堆積角試驗過程如圖5所示，測得實際堆積角平均值為31.67°。

圖5 堆積角左邊緣擬合 Fig.5 Left edge fitting of actual repose angle

2.2 油菜莖稈模型建立及EDEM參數設置

在EDEM中采用Hertz-Mindlin無滑移模型進行油菜莖稈堆積試驗（圖6），為簡化模型并提高仿真效率，采用圓形顆粒組合的方式形成與物理試驗相同長度和直徑的油菜莖稈，直徑10 mm莖稈顆粒模型如圖6a所示。通過多次預仿真試驗并結合文獻[14，34]及前文試驗適當擴大試驗參數取值范圍，如表3所示。

表3 油菜莖稈顆粒堆積角仿真模型參數 Table 3 Parameters of repose angle simulation model for rape stalk particles

在EDEM仿真試驗中，建立與物理試驗參數相同的鋼制圓筒與底板模型，在圓筒上方建立Polygon平面作為顆粒工廠，經過多次重復試驗，設定每種顆粒生成速率為75 粒/s，共生產時間為2 s，顆粒產生完全并靜置穩定顆粒堆后采用與物理試驗相同的速度進行提升；仿真總時間為4.5 s，時間步長為Rayleigh時間步長的20%，網格尺寸為顆粒半徑的3倍，EDEM仿真試驗結果如圖6b所示，仿真試驗圖像處理如圖6c所示。

圖6 堆積角仿真試驗 Fig.6 Simulation test of repose angle

2.3 仿真試驗

為減少試驗次數，應用Design Expert 10.0.7進行Plackett-Burman篩選試驗，以莖稈堆積角為響應值，篩選出對堆積角影響顯著的參數。將表3中的～進行最大、最小值分別編碼水平+1、-1，編碼結果如表4所示。分別測量左右兩側的堆積角數值，然后取平均值，Plackett-Burman試驗結果如表5所示。

表4 Plackett-Burman試驗方案 Table 4 Plackett-Burman test scheme

表5 Plackett-Burman試驗結果 Table 5 Plackett-Burman test results

用Design Expert 10.0.7進行方差分析，得到各個參數的影響效果，如表6所示。根據貢獻率進行顯著性排序，在油菜莖稈堆積試驗中，，對堆積角的影響較為顯著，因此只考慮這3個因素進行最陡爬坡試驗；其余顯著性影響較小的取值分別為：取0.4，取1.5 MPa，取0.277，取0.790，取0.391。

表6 Plackett-Burman試驗參數顯著性分析 Table 6 Significance analysis of Plackett-Burman test parameters

基于Plackett-Burman堆積試驗篩選出3個顯著性影響參數（莖稈-莖稈靜摩擦系數，莖稈-莖稈動摩擦系數，莖稈-鋼板動摩擦系數），以仿真堆積角與實際堆積角的相對誤差作為評價指標來確定試驗參數的最佳范圍，試驗方案及結果如表7所示。

表7 最陡爬坡試驗設計方案及結果 Table 7 The steepest climbing test design scheme and results

3個試驗因素逐漸遞增，堆積角相對誤差呈先減小后增加的變化趨勢，其中4號的堆積角相對誤差最小，為1.89%，基于最陡爬坡試驗確定以4號試驗中的各個參數作為后期試驗的中心點，3號、5號作為低水平與高水平，通過Box-Behnken響應面試驗驗證該模型的顯著性關系，并通過優化得出最佳參數組合。由式（7）計算出仿真堆積角與實際堆積角的相對誤差。

利用Design-Expert 10.0.7進行三因素三水平響應曲面試驗設計，共進行17組仿真試驗，試驗設計方案與結果如表8所示。對Box-Behnken試驗結果進行多元回歸擬合，得到油菜莖稈堆積角與3個顯著性參數二階回歸方程為

表8 Box-Behnken試驗設計方案及結果 Table 8 Box-Behnken experiment design scheme and results

由該模型方差分析結果可知（表9所示），莖稈-莖稈動摩擦系數對堆積角影響極顯著；莖稈-莖稈靜摩擦系數，莖稈-鋼板動摩擦系數對堆積角影響顯著；、、對堆積角影響極顯著；對堆積角影響顯著；、對堆積角影響不顯著；該模型的＜0.05表明該模型的各個參數與響應值之間的關系顯著；失擬項=0.356 2＞0.05，表明該模型擬合良好；變異系數（CV）為1.62%，說明本試驗具有較好的可靠性；決定系數=0.966 6，校正決定系數=0.923 8，表明擬合方程可靠度較高；精密度（Adeq Precisior）為15.726，表明該模型具有良好的精確度。

表9 Box-Behnken試驗回歸模型方差分析 Table 9 Variation analysis of Box-Behnken quadratic model

根據回歸模型方差分析可知，、這兩個交互項對堆積角的影響不顯著，交互項對堆積角的影響極顯著。應用Design-Expert 10.0.7軟件繪制三個顯著性因素交互作用堆積角響應曲面，如圖7所示。由圖7a可知，當莖稈-莖稈間靜摩擦系數一定時，隨著莖稈-莖稈動摩擦系數的增加，堆積角變化趨勢明顯，對堆積角影響極顯著，而莖稈-莖稈動摩擦系數一定時，隨著莖稈-莖稈間靜摩擦系數的增加，堆積角呈先增后減的趨勢，對堆積角影響較小；由圖7b可知，響應曲面為凸型曲面，莖稈-莖稈靜摩擦系數和莖稈與鋼板動摩擦系數對堆積角的影響均呈先增后減趨勢，表明兩因素對堆積角的影響基本相似；由圖7c可知，莖稈與莖稈動摩擦系數對響應面趨勢較陡，該因素比莖稈與鋼板動摩擦系數影響顯著。

圖7 物料參數之間的交互作用響應曲面 Fig.7 Response surfaces of interaction between material parameters

利用Design-Expert 10.0.7軟件的優化模塊中以實際堆積角為目標值（31.67°）進行求解，對回歸模型進行尋優，得到一組與物理試驗測定結果較為相近的一組解為：莖稈-莖稈靜摩擦系數為0.707，莖稈-莖稈動摩擦系數為0.015，莖稈-鋼板動摩擦系數為0.012，其他非顯著性參數取平均值。為驗證離散元仿真的可靠性，以上述參數作為EDEM仿真參數，進行3組仿真試驗，得到的休止角為30.90°、31.81°、31.80°，平均值為31.50°，與物理試驗結果的相對誤差為0.54%，驗證了仿真試驗的真實性與可靠性。

3 旋風分離清選仿真試驗及驗證

3.1 DEM-CFD標定參數驗證

利用 Solidworks對旋風分離筒進行建模，在Workbench 17.0的mesh中采用四面體非結構化方法劃分旋風分離器的網格，對進出口處網格加密處理，并進行邊界層網格的劃分，以保證其計算結果準確性，總體網格數約75萬（如圖8a）。對莢殼的模型進行三維掃描后在EDEM中進行實心顆粒填充；油菜籽粒簡化為2 mm的球體；測量凹板篩分離出的油菜脫出物中短莖稈長度與直徑，為簡化計算量，結合文獻[34-36]對莖稈顆粒進行簡化，采用較為普遍應用的長度40 mm、直徑3 mm的圓柱體短莖稈，如圖8b所示。

圖8 仿真網格及顆粒模型 Fig.8 Simulation grid and particle models

進料口為速度入口，吸雜口、出料口均為壓力出口。EDEM和CFD的時間步長分別設置為5×10和1×10，氣流介質為空氣，其密度和黏性系數分別為1.225 kg/m和1.7894×10kg/(m·s)，重力加速度為9.81 m/s。物料顆粒力學參數及接觸參數（非顯著性參數取試驗平均值）如表10和表11所示。旋風分離筒涉及復雜的密相氣體-顆粒兩相流動，在進行DEM-CFD耦合試驗時，采用考慮顆粒對流體影響的Eulerian模型；顆粒相由EDEM 2018求解，為簡化計算量、提高仿真效率，采用Hertz-Mindlin 無滑移接觸模型，可以較好反應剛性物體間的力學行為；氣相采用Fluent 17.0求解，Fluent將某一時間的流場計算至收斂后，通過曳力模型將流場信息轉化為EDEM中作用在顆粒上的流體曳力，EDEM計算此時每個顆粒所受的外力（流體曳力、碰撞力等）從而更新顆粒此時的位置、速度、加速度等信息，最后這些信息以動量匯形式加載到CFD中，實現相間的耦合，由于清選筒內是三維強旋流流場，故選擇RNG-模型，勾選適用于旋流的Swirl Dominated Flow選項。

表10 物料力學參數 Table 10 Mechanical parameters of materials

表11 物料接觸參數 Table 11 Contact coefficient of materials

旋風分離筒內的氣流分布規律及物料間的接觸參數決定了筒內物料的運動狀態，調節旋風分離筒的透明度，完成可視化的清選分離仿真過程顯示。由圖9可以看出物料在筒內的運動狀態，物料從切向進入旋風分離筒內并繞著筒內壁向下運動，莖稈和莢殼在運動過程中被吸雜口處的負壓吸走，籽粒沿著內壁不斷向下運動從出料口處流出或者部分積聚在出料口的底部，總體來看數值模擬可以較為準確的反應籽粒及雜余的運動軌跡。在EDEM后處理模塊Grid Bin Group統計出料口和吸雜口處的籽粒、莢殼、短莖桿的質量，通過式（9）計算得到旋風分離清選裝置油菜籽粒清潔率為94.42%，損失率為3.96%。

圖9 筒內物料運動狀態及運動軌跡 Fig.9 Movement state and motion trajectory of the material in the cylinder

式中Y為油菜籽粒清潔率，%；為出料口油菜籽粒質量，g；為出料口油菜脫出物總質量，g；Y為損失率，%；為吸雜口油菜籽粒質量，g。

3.2 臺架驗證試驗

本試驗依托華中農業大學工學院自主研發的收獲關鍵部件試驗臺開展臺架驗證試驗，試驗臺主要包括輸送帶、割臺、縱軸流脫粒分離裝置和旋風分離清選裝置等組成，如圖10所示。收獲關鍵部件試驗臺相關參數實時可調并且相關參數可以在線監測。

圖10 臺架試驗 Fig.10 Bench test

試驗結束后，收集出料口及吸雜口處物料進行人工分選并稱量，通過式（9）計算油菜籽粒清潔率與損失率。入料口速度為6 m/s，吸雜口速度為13 m/s時，旋風分離清選裝置油菜籽粒清潔率為91.84%，損失率為4.28%。仿真試驗與臺架試驗的油菜籽粒清潔率與損失率相對誤差分別為2.81%和7.48%。

4 結 論

1）研究得出了油菜聯合收獲脫出物的相關基礎參數。開展了油菜脫出物特征參數、接觸參數測量，測得油菜莖稈的泊松比、剪切模量分別為0.4、1.5 MPa；鋼板與籽粒、莖稈、莢殼的靜摩擦系數與碰撞恢復系數范圍分別為0.214～0.390、0.711～0.869、0.647～0.781；0.543～0.667、0.212～0.570、0.141～0.283。莖稈與莖稈、莢殼的靜摩擦系數與碰撞恢復系數范圍分別為0.650～0.754、0.637～0.751；0.187～0.367、0.158～0.265。莖稈與籽粒的靜摩擦系數與碰撞恢復系數分別為0.292～0.516；0.515～0.632；籽粒與莢殼的靜摩擦系數與碰撞恢復系數范圍分別為0.635～0.907；0.566～0.612，莢殼與莢殼的靜摩擦系數與碰撞恢復系數范圍分別為0.625～0.724；0.141～0.265。

2）開展了油聯合收獲油菜莖稈離散元參數標定，通過Plackett-Burman試驗，篩選出對堆積角影響顯著的參數為：莖稈-莖稈靜摩擦系數、莖稈-莖稈動摩擦系數、莖稈-鋼板動摩擦系數，通過Design-Expert 10.0.7軟件優化模塊得到顯著性參數最佳組合參數為0.707、0.015、0.012，在最優組合參數下仿真堆積角平均值為31.50°，與實際堆積角的相對誤差為0.54%。

3）結合標定的油菜脫出物不同組分間接觸參數，開展了基于DEM-CFD旋風分離清選過程數值分析和臺架驗證試驗。結果表明：旋風分離清選仿真試驗油菜籽粒清潔率為94.42%，損失率為3.96%；仿真試驗與臺架試驗籽粒清潔率與損失率相對誤差分別為2.81%和7.48%。表明標定的油菜脫出物接觸參數可為油菜聯合收獲機具的離散元仿真分析提供基礎參數。

聯合收獲油菜莖稈取材和參數標定還需進一步擴大樣本，后續將持續開展不同地區、不同部位油菜莖稈力學特性的具體性研究；莢殼易碎、呈長扁內凹壁結構，常規力學測試方法難以測出其特征參數和接觸參數；后續將結合高速攝影技術和能量守恒定律對莢殼及脫出物不同組分間接觸參數進一步標定。