南夏豆25收獲期植株離散元建模方法與參數標定

譚云峰，陳 霖，胡 森，王 鍵，陳治帆，呂小榮

（四川農業大學機電學院，四川 雅安 625014）

大豆作為我國重要農業作物之一，是植物油和植物蛋白主要來源，對我國油脂安全供給、畜牧業發展及國家糧食安全具有重要意義[1-2]。大豆機械化收獲是實現大豆生產全程機械化重要環節。作為收獲機核心部件，脫粒裝置對整機工作質量起決定性作用[3]，收獲過程損失降低經濟效益。莖稈是大豆收獲過程中主要雜質，其物理特性對收獲過程影響較大。脫粒過程復雜，試驗中較難獲得物料狀態變化數據。因此，有必要建立還原度較高的大豆莖稈模型，為脫粒過程中莖稈與脫粒裝置交互作用研究提供可靠依據。

目前已有針對植株仿真試驗的相關研究。大多采用圓形顆粒堆砌成剛性長桿模型，但植株脫粒過程中出現彎曲變形、破碎等現象，傳統剛性植株模型在仿真計算過程中被視作無任何位移變形的剛體。劉基、謝干等使用莖稈幾何參數建立剛性大豆、小麥莖稈模型，利用離散元仿真研究脫粒過程中籽粒破碎情況、功耗影響因素[4-5]。王萬章等在多球面填充法構建剛性莖稈模型基礎上，基于Hertz－Mindlin with bonding接觸模型構建小麥柔性莖稈模型，該接觸模型可承受一定量法向和切向位移，直至到達臨界法向和切向力，模型可彎曲和破碎[6]；趙吉坤等利用EDEM 軟件中Hertz-Mindlin with Bonding 模型建立水稻秸稈模型，在秸稈彎曲力學特性試驗基礎上，對模型參數進行標定[7]；Wang 等利用空心圓柱彈性鍵建立水稻植株離散元模型并在顆粒水平上準確描述柔性水稻植株在多重相互作用條件下動態行為[8]；Mao 等使用空心圓柱鍵建立柔性管模型并驗證可有效再現模型流動行為，用于建立新秸稈模型[9]；Leblicq 等建立DEM中莖稈彎曲模型，可用于單個莖稈彎曲行為驗證[10]。Yu 等基于離散元方法，建立玉米穗數值模型并用于玉米脫粒過程仿真，仿真結果與試驗數據吻合[11-12]。張鋒偉等早期研究結果表明，分析玉米秸稈破碎揉絲過程中受力，可表征玉米秸稈力學特性，結合力學特性試驗對玉米秸稈進行顆粒黏結模型黏結參數校核[13]。劉禹辰等通過對玉米秸稈進行壓縮、玉米秸稈外表皮與內穰拉伸和剪切等力學試驗，采用離散元法進行力學特性仿真，將黏結鍵斷裂情況作為玉米秸稈受外力時破裂程度評定指標[14]。張開飛等針對大豆秸稈進行彎曲、剪切力學試驗發現，載荷加載速度對大豆秸稈彎曲應力及剪切應力產生影響，同一株大豆秸稈不同部位其表現強度也有較大差異[15]。?ztürk等研究大豆莖稈切削性能，切削力隨莖稈直徑增加而線性增加，不同高度間莖稈切削力無顯著差異[16]?；陔x散元法建立的莖稈模型已廣泛應用于現代農業裝備研究，但大豆莖稈柔性離散元模型建立與黏結參數標定研究仍未見報道。

針對實際脫粒過程中大豆莖稈出現彎曲、破碎等情況，本文以南夏豆25（其在四川省平壩、丘陵及低山區廣泛種植）為研究對象，通過彎曲物理試驗測得大豆莖稈參數極限斷裂載荷和極限斷裂位移，利用EDEM 軟件中Hertz-Mindlin with bond?ing 接觸模型，構建可破碎柔性大豆植株離散元模型，以極限斷裂載荷和極限斷裂位移為響應值，將其彎曲仿真試驗與物理試驗對比，分析大豆Bonding 參數及其交互作用，并對其Bonding 參數進行標定。以籽粒損失率為指標，進行仿真與物理對比試驗，驗證模型準確性。

1 材料與方法

1.1 大豆植株基本物性參數確定

本研究大豆植株取樣于四川省仁壽縣玉米-大豆間套作種植基地，品種為南夏豆25，為獲取大豆植株幾何尺寸，利用游標卡尺（精度：0.01 mm）進行測量，樣品數量為100株，結果如表1所示。

表1 大豆植株幾何尺寸Table 1 Plant geometry of soybean

1.2 大豆莖稈物理彎曲試驗

收獲期大豆莖稈含水率為40%~50%，但大豆莖稈含水率對大豆莖稈黏聚力、流動函數值并無顯著影響，僅對其自身力學物理性質有影響[17]。為測得大豆莖稈極限破碎位移和極限破碎載荷，選用直徑為7.13~7.75 mm莖稈為試驗材料，在三點彎曲夾具質構儀（精度為0.001 mm、0.001 N、刀刃直徑為10 mm，購自上海保圣實業發展有限公司）上作彎曲試驗，試驗過程見圖1。試驗重復20次取平均值，得到大豆莖稈極限位移和載荷，位移載荷曲線見圖2。由圖2可見，隨刀刃下壓，莖稈開始彎曲發生彈性變形，大豆莖稈所承受載荷迅速增大，當達到最大載荷98.994 N 時，大豆莖稈加載位移為3.559 mm，莖稈與儀器接觸點四周產生縱向裂橫，伴隨斷面增大的同時，載荷力減小，直到莖稈下表面出現橫向斷裂，載荷下降至零。

圖1 大豆植株彎曲過程Fig.1 Bending process of soybean plant

圖2 大豆植株彎曲力-位移曲線變化Fig.2 Change of curve of bending force-displacement of soybean plant

2 收獲期大豆植株離散元仿真模型創建與仿真試驗

2.1 離散單元法基本原理及接觸模型選取

離散單元法最初由美國學者Cundall P.A 在20世紀70 年代初提出。離散單元法是將不連續體分離為剛性單元體，各單元體滿足運動方程，運用迭代方法進行循環迭代計算剛性單元體的運動方程，得到研究對象宏觀運動規律[18-21]。

Hertz-Mindlin with bonding V2接觸模型采用球體之間鑲嵌重疊將顆粒組合體連接，球體單元以Bond 鍵相互粘結方式構造，黏結點可承受來自外界的切向和法向位移，黏結斷裂需滿足達到最大法向和切向剪切應力。Hertz-Mindlin with bonding V2 接觸模型適用于顆粒破碎、分解的問題。根據大豆莖稈可彎曲、破碎特性，采用該接觸模型研究大豆植株離散元參數，其模型原理圖和Bond 鍵示意分別如圖3和4所示。

圖3 Bonding 模型原理圖Fig.3 Schematic diagram of the Bonding model

圖4 Bond鍵示意圖Fig.4 Bond bond diagram

2.2 大豆植株黏結離散元模型建立

通過前期預試驗，得出收獲期不同含水率豆莢破莢力分布在0.76~6.64 N，脫粒作業中脫粒齒提供打擊力為40~70 N，遠高于豆莢破莢力，實際脫粒作業中豆莢破開籽粒離開豆莢即視為完成脫粒作業，為簡化模型，可將豆莢簡化為籽粒?？紤]仿真模型計算體量，將大豆植株視為圓柱體，基于Hertz-Mindlin with bonding V2 黏結接觸模型，利用半徑3.5 mm 球形粒子相互黏結構建大豆莖稈離散元Bonding V2模型。大豆籽粒形狀與橢球體相似，因此在建立顆粒模型時將大豆種子簡化為橢球體，測得大豆籽粒尺寸數據后，利用五球填充方法構建大豆籽粒仿真模型，可反映大豆籽粒接觸狀態[22-23]。為分析脫粒裝置脫粒分離能力和夾帶損失率，籽粒與莖稈間同樣可利用Bond 鍵黏結。籽粒與莖稈質量比例等于喂入物料草谷比3∶2，收獲時期大豆植株離散元模型如圖5所示。

圖5 收獲期大豆植株離散元模型Fig.5 Discrete element model of soybean plant at harvest

2.3 仿真試驗參數設置

利用離散元仿真軟件EDEM 2020 建立與物理試驗一致的大豆莖稈彎曲數值模擬試驗，如圖6所示。仿真試驗數據保存間隔為0.001 s，總時長預設為10 s，時間步長為8.7846×10-7s，網格最小顆粒半徑劃分為3倍。通過預試驗及文獻得到仿真中用到的材料及接觸參數[24-25]，數值如表2所示。

圖6 大豆莖稈彎曲數值模擬試驗Fig.6 Numerical simulation experiment of soybean stem bending

表2 離散元仿真參數Table 2 Discrete element simulation parameters

2.4 試驗方?法設計

Bonding 模型由單位面積法向剛度、切向剛度、臨界法向應力、臨界切向應力和黏結半徑比例5個參數決定最終黏結鍵強度。為探究各參數對極限位移及載荷影響的顯著性，以極限斷裂位移Y1和極限斷裂載荷Y2為目標，以單位面積法向剛度、單位面積切向剛度、臨界法向應力、臨界切向應力和黏結半徑比例為標定參數。為獲取對試驗結果有顯著效應的參數，采用Plackett-Burman試驗方法，觀察顯著因素t值正負效應，確定最陡爬坡試驗路徑方向，對仿真過程中其余不顯著因素取中間值；基于上述Plackett-burman 試驗結果，進行最陡爬坡試驗，確定顯著仿真因素最佳范圍；最后采用Box-Behnken試驗進行響應面優化分析，獲取試驗仿真最佳參數。

3 仿真試驗結果與分析

3.1 Plackett-burman試驗

Plackett-burman 試驗5 個因子為單位面積法向剛度、單位面積切向剛度、臨界法向應力、臨界切向應力和黏結半徑（X1、X2、X3、X4、X5），選用N=11 的Plackett-burman 試驗設計，預留6 個虛擬項作誤差分析，試驗因素水平如表3所示。

表3 Plackett-burman試驗因素水平Table 3 Plackett-burman trial factor levels

Plackett-burman 試驗結果如表4所示，為得到各因素顯著效應情況，利用Design-expert 11.0 進行方差分析和t檢驗選取影響顯著的因素，結果分別如表5、圖7所示。

圖7 帕累托圖Fig.7 Pareto chart

表4 Plackett-burman試驗結果Table 4 Plackett-burman trial results

表5 Plackett-burman試驗結果分析Table 5 Plackett-burman test results analysis

通過表5可知，對極限斷裂位移Y1的影響因素從大到小依次為X3、X1、X5、X2、X4，其中顯著影響因素為X1、X3，由圖7 帕累托圖（Pareto Chart）的t-value檢驗不僅可得到各因素顯著性排序，還可觀察到因素正負效應，其中X1、X5、X2對目標值影響為負效應，X3、X4為正效應；對極限斷裂載荷Y2的影響因素從大到小依次為X5、X3、X1、X4、X2，其中影響顯著因素為X1、X3、X5，且X1、X4對目標值的影響為負效應，X3、X5對目標值的影響為正效應。

剔除對目標值影響不顯著因素，各因素與極限斷裂位移Y1和極限斷裂載荷Y2的模型為：

3.2 最陡爬坡試驗

根據Plackett-burman試驗結果可知，X2、X4對目標值影響均不顯著，因此X2值為1.44×109N·m-3、X4值為1.26×106Pa。選取因素X1、X3、X5初值分別為3.6×109N·m-3、2.1×106Pa、1.6，步長分別為-1.5×108N·m-3、1.8×105Pa、0.16 進行最陡爬坡試驗，以求出更靠近標準值的參數組合，最陡爬坡試驗方案設計及結果如表6所示。由前文可知，莖稈彎曲破壞時極限斷裂位移和極限斷裂載荷分別為3.559 mm、98.944 N，表6 仿真試驗結果中位移及載荷與響應值的誤差均先減后增。其中第4組試驗結果與響應值誤差最小，分別為0.67%、3.63%，因此最佳參數組合取值范圍在試驗4附近。繼續取試驗3、試驗4、試驗5 參數組合分別為低、中、高水平進行后續Box-Behnken數值模擬試驗。

表6 最陡爬坡試驗方案及結果Table 6 Program and results of the steepest climbing test

3.3 Box-Behnken試驗

Plackett-burman試驗篩選出具有顯著效應的參數，最陡爬坡試驗確定最佳范圍，現以試驗3、試驗4、試驗5參數組合分別為低（-1）、中（0）、高水平（1）進行Box-Behnken 試驗，得出試驗因素水平設計，如表7所示。試驗選取各中心點對試驗誤差進行評估，試驗方案及結果如表8所示。方差分析結果如表9所示。

表7 試驗因素水平Table 7 Level of test factors

表8 Box-Behnken試驗方案及結果Table 8 Box-behnken test protocol and results

表9 Box-Behnken試驗方差分析Table 9 Box-behnken test analysis of variance

由表9可知，3個因素對極限斷裂位移Y1和極限斷裂載荷Y2的二次回歸模型均顯著且其失擬項不顯著，表明能夠準確預測兩個目標值變化趨勢。其中對Y1影響顯著性由大到小為X3、X5、X1X5、X3、X1，其中X1不顯著。

為得到更為精確的參數組合，需縮小3個因素取值范圍。已知Y1和Y2目標值分別為3.559 mm、98.944 N，X2取值為1.44×109N·m-3，X4取值為1.26×106Pa，利用Design-expert 11.0 軟件對試驗因素進行優化[26-27]。得到最優參數組合X1、X3、X5分別為3.1277878×109N·m-3、2.657545×106Pa、2.1，此時Y1、Y2分別為3.575 mm、99.12 N，與目標值相對誤差分別為0.45%、0.13%。

4 參數優化試驗驗證

4.1 試驗條件

為進一步驗證優化得到的大豆植株Bonding 模型黏結參數組合正確性，選取4LZ-1.6Z 型大豆收獲機上常用縱軸流脫粒裝置進行脫粒性能試驗，搭建試驗臺架如圖8a 所示，試驗材料為提前在田間采收的南夏豆25 植株。每組試驗設定喂入量為1.5 kg·s-1，總喂入時間2 s，在不同脫粒轉速下試驗，每組試驗重復3 次取平均值。采用3D 建模軟件Solidworks 對脫粒裝置進行建模，如圖8b 所示。試驗選取脫出物夾帶損失率為試驗指標，脫出物收集盤位置如圖9所示，具體計算方法如下：式中，Y3-夾帶損失率（%）；m1-前收集盤中收集籽?？傎|量（g）；m2-后收集盤中收集籽?？傎|量（g）。

圖8 脫粒臺架試驗Fig.8 Threshing bench test

圖9 收集盤位置示意Fig.9 Collection disk position

4.2 試驗結果

在滾筒轉速分別為350、400、450 r·min-1時得到脫粒裝置物理試驗與仿真試驗夾帶損失率如圖10 所示?？梢?，仿真試驗夾帶損失率與物理試驗平均誤差為7.07%，表明本文標定的Bonding 模型參數組合準確度高。

圖10 仿真試驗與臺架試驗得損失率對比Fig.10 Comparison of loss rate between simulation test and bench test

籽粒分布情況如圖11所示。為更好統計大豆籽粒軸向分布規律，對前收集盤作軸向編號如圖12所示，籽粒軸向分布數據如圖13 所示?？梢?，大豆籽粒經脫粒裝置脫粒分布大致相同，均呈遞減趨勢。

圖11 籽粒分布情況Fig.11 Grain distribution

圖12 前收集盤軸向編號Fig.12 Front collection disc axial number

圖13 籽粒質量軸向分布Fig.13 Axial distribution of grain mass

大豆植株離散元模型在脫粒元件與凹板篩打擊、搓擦作用下，跟隨釘齒旋轉同時進行軸向流動，當外加載荷超過黏結鍵強度時產生黏結鍵斷裂，物料運動過程如圖14 所示，表示物料在不同時刻狀態，包括在凹板篩低處壓縮及導流板處拋散。

圖14 脫粒過程中大豆莖稈運動過程Fig.14 Movement of soybean stem during threshing

由圖14 可見，物料隨脫粒元件流動同時，易堆積在凹板篩右下部，此時物料受到打擊、搓擦作用最強，最易破碎，造成物料破碎率高、脫粒功耗增大，導致含雜率增加、脫粒滾筒易堵塞。仿真情況與物理試驗時情況一致，進一步驗證本文標定的仿真參數準確性。

5 結 論

a.為測得所脫物料力學性能，利用質構儀與三點彎曲試驗探頭，通過彎曲試驗得到大豆植株極限斷裂位移和極限斷裂載荷分別為3.559 mm、98.944 N。

b. 為對所建植株模型粘結參數進行標定，利用Placket-burman試驗找出顯著因子、最陡爬坡試驗縮小參數區間、Box-Behnken試驗確定最優參數組合。Placket-burman 試驗結果表明，對大豆莖稈極限斷裂位移有顯著影響的因素為：單位面積法向剛度和臨界法向應力；對極限斷裂載荷有顯著的影響因素為：單位面積法向剛度、臨界法向應力和黏結半徑比例。Box-Behnken 試驗結果表明，單位面積法向剛度、臨界法向應力和粘結半徑比例對極限斷裂位移影響顯著，臨界法向應力和粘結半徑比例對極限斷裂載荷影響顯著。通過優化工具求解得到當單位面積法向剛度、切向剛度、臨界法向應力、切向應力、粘結半徑比例分別為3.1277878×109N·m-3、1.44×109N·m-3、2.657545×106Pa、1.26×106Pa、2.1時，極限斷裂位移和極限斷裂載荷分別為3.575 mm、99.12 N，與真實值相對誤差分別為0.45%、0.13%。

c.為驗證植株模型準確性，以相同條件進行大豆脫粒臺架試驗，對比仿真和物理試驗損失率，兩者平均誤差為7.07%，證明所提出大豆植株Bonding模型參數可靠。