南夏豆25收獲期植株離散元建模方法與參數(shù)標(biāo)定

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

（四川農(nóng)業(yè)大學(xué)機(jī)電學(xué)院，四川 雅安 625014）

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

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

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

1 材料與方法

1.1 大豆植株基本物性參數(shù)確定

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

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

1.2 大豆莖稈物理彎曲試驗(yàn)

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

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

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

2 收獲期大豆植株離散元仿真模型創(chuàng)建與仿真試驗(yàn)

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

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

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

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

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

2.2 大豆植株黏結(jié)離散元模型建立

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

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

2.3 仿真試驗(yàn)參數(shù)設(shè)置

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

圖6 大豆莖稈彎曲數(shù)值模擬試驗(yàn)Fig.6 Numerical simulation experiment of soybean stem bending

表2 離散元仿真參數(shù)Table 2 Discrete element simulation parameters

2.4 試驗(yàn)方?法設(shè)計(jì)

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

3 仿真試驗(yàn)結(jié)果與分析

3.1 Plackett-burman試驗(yàn)

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

表3 Plackett-burman試驗(yàn)因素水平Table 3 Plackett-burman trial factor levels

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

圖7 帕累托圖Fig.7 Pareto chart

表4 Plackett-burman試驗(yàn)結(jié)果Table 4 Plackett-burman trial results

表5 Plackett-burman試驗(yàn)結(jié)果分析Table 5 Plackett-burman test results analysis

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

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

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

根據(jù)Plackett-burman試驗(yàn)結(jié)果可知，X2、X4對目標(biāo)值影響均不顯著，因此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 進(jìn)行最陡爬坡試驗(yàn)，以求出更靠近標(biāo)準(zhǔn)值的參數(shù)組合，最陡爬坡試驗(yàn)方案設(shè)計(jì)及結(jié)果如表6所示。由前文可知，莖稈彎曲破壞時(shí)極限斷裂位移和極限斷裂載荷分別為3.559 mm、98.944 N，表6 仿真試驗(yàn)結(jié)果中位移及載荷與響應(yīng)值的誤差均先減后增。其中第4組試驗(yàn)結(jié)果與響應(yīng)值誤差最小，分別為0.67%、3.63%，因此最佳參數(shù)組合取值范圍在試驗(yàn)4附近。繼續(xù)取試驗(yàn)3、試驗(yàn)4、試驗(yàn)5 參數(shù)組合分別為低、中、高水平進(jìn)行后續(xù)Box-Behnken數(shù)值模擬試驗(yàn)。

表6 最陡爬坡試驗(yàn)方案及結(jié)果Table 6 Program and results of the steepest climbing test

3.3 Box-Behnken試驗(yàn)

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

表7 試驗(yàn)因素水平Table 7 Level of test factors

表8 Box-Behnken試驗(yàn)方案及結(jié)果Table 8 Box-behnken test protocol and results

表9 Box-Behnken試驗(yàn)方差分析Table 9 Box-behnken test analysis of variance

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

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

4 參數(shù)優(yōu)化試驗(yàn)驗(yàn)證

4.1 試驗(yàn)條件

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

圖8 脫粒臺(tái)架試驗(yàn)Fig.8 Threshing bench test

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

4.2 試驗(yàn)結(jié)果

在滾筒轉(zhuǎn)速分別為350、400、450 r·min-1時(shí)得到脫粒裝置物理試驗(yàn)與仿真試驗(yàn)夾帶損失率如圖10 所示?？梢姡抡嬖囼?yàn)夾帶損失率與物理試驗(yàn)平均誤差為7.07%，表明本文標(biāo)定的Bonding 模型參數(shù)組合準(zhǔn)確度高。

圖10 仿真試驗(yàn)與臺(tái)架試驗(yàn)得損失率對比Fig.10 Comparison of loss rate between simulation test and bench test

籽粒分布情況如圖11所示。為更好統(tǒng)計(jì)大豆籽粒軸向分布規(guī)律，對前收集盤作軸向編號(hào)如圖12所示，籽粒軸向分布數(shù)據(jù)如圖13 所示?？梢姡蠖棺蚜＝?jīng)脫粒裝置脫粒分布大致相同，均呈遞減趨勢。

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

圖12 前收集盤軸向編號(hào)Fig.12 Front collection disc axial number

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

大豆植株離散元模型在脫粒元件與凹板篩打擊、搓擦作用下，跟隨釘齒旋轉(zhuǎn)同時(shí)進(jìn)行軸向流動(dòng)，當(dāng)外加載荷超過黏結(jié)鍵強(qiáng)度時(shí)產(chǎn)生黏結(jié)鍵斷裂，物料運(yùn)動(dòng)過程如圖14 所示，表示物料在不同時(shí)刻狀態(tài)，包括在凹板篩低處壓縮及導(dǎo)流板處拋散。

圖14 脫粒過程中大豆莖稈運(yùn)動(dòng)過程Fig.14 Movement of soybean stem during threshing

由圖14 可見，物料隨脫粒元件流動(dòng)同時(shí)，易堆積在凹板篩右下部，此時(shí)物料受到打擊、搓擦作用最強(qiáng)，最易破碎，造成物料破碎率高、脫粒功耗增大，導(dǎo)致含雜率增加、脫粒滾筒易堵塞。仿真情況與物理試驗(yàn)時(shí)情況一致，進(jìn)一步驗(yàn)證本文標(biāo)定的仿真參數(shù)準(zhǔn)確性。

5 結(jié) 論

a.為測得所脫物料力學(xué)性能，利用質(zhì)構(gòu)儀與三點(diǎn)彎曲試驗(yàn)探頭，通過彎曲試驗(yàn)得到大豆植株極限斷裂位移和極限斷裂載荷分別為3.559 mm、98.944 N。

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

c.為驗(yàn)證植株模型準(zhǔn)確性，以相同條件進(jìn)行大豆脫粒臺(tái)架試驗(yàn)，對比仿真和物理試驗(yàn)損失率，兩者平均誤差為7.07%，證明所提出大豆植株Bonding模型參數(shù)可靠。