水稻莖稈接觸物理參數(shù)測定與離散元仿真標(biāo)定

侯 杰，謝方平，2*，王修善，2，劉大為，2，陳志剛

（1.湖南農(nóng)業(yè)大學(xué) 機(jī)電工程學(xué)院，湖南 長沙 410128；2.湖南省智能農(nóng)機(jī)裝備重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410128）

【研究意義】水稻是我國主要糧食作物之一，其種植面積廣、產(chǎn)量大[1]，機(jī)械化收獲水平逐年提高[2]，但是由于水稻生物學(xué)特性變化，水稻機(jī)械化收獲清選過程中物料堵塞問題嚴(yán)重[3]，解決該問題的前提是建立較為精確的水稻模型，模擬其在聯(lián)合收割機(jī)中的運(yùn)動(dòng)規(guī)律，深入研究堵塞原因，可為后續(xù)聯(lián)合收割機(jī)優(yōu)化設(shè)計(jì)提供理論基礎(chǔ)。【前人研究進(jìn)展】離散元法是專門用來解決不連續(xù)介質(zhì)問題的數(shù)值模擬方法[3]，在離散元仿真模擬中，參數(shù)標(biāo)定是最關(guān)鍵的問題之一，EDEM 模擬谷物脫粒與清選時(shí)[4-9]，莖稈等的參數(shù)設(shè)置直接影響仿真結(jié)果的可信度。近年來，離散元參數(shù)標(biāo)定方法不斷增多[10]，劉凡一等[11]采用三點(diǎn)彎曲試驗(yàn)對小麥短莖稈粘結(jié)參數(shù)進(jìn)行離散元標(biāo)定，結(jié)果表明：莖稈的彈性模量隨著粘結(jié)半徑、單位面積法向剛度和單位面積剪切剛度的增大而增大；張濤等[12]以徑向堆積角相對誤差值為評價(jià)指標(biāo)，應(yīng)用正交方法標(biāo)定玉米秸稈離散元仿真需要輸入的接觸參數(shù)；王韋韋等[13]通過堆積角試驗(yàn)和?？讐嚎s試驗(yàn)標(biāo)定玉米秸稈接觸參數(shù)，試驗(yàn)結(jié)果表明：仿真試驗(yàn)堆積角與物理堆積角相對誤差為0.68%，仿真試驗(yàn)與物理試驗(yàn)?zāi)？讐嚎s比相對誤差為0.98%，表明標(biāo)定的離散元參數(shù)可用作玉米秸稈離散元仿真研究；馬彥華等[14]以休止角為目標(biāo)值對苜蓿秸稈壓縮過程中離散元仿真研究所用接觸參數(shù)進(jìn)行標(biāo)定；王萬章等[15]通過斜面法對比仿真與物理試驗(yàn)，對小麥莖稈接觸參數(shù)進(jìn)行標(biāo)定；廖宜濤等[16]通過堆積角和響應(yīng)面試驗(yàn)標(biāo)定了油菜莖稈離散元接觸和粘結(jié)參數(shù)。

【本研究切入點(diǎn)】較多研究者對莖稈離散元參數(shù)標(biāo)定提供了可行的方法，但少有對水稻莖稈接觸參數(shù)和粘結(jié)參數(shù)的全部標(biāo)定，且試驗(yàn)設(shè)計(jì)所用方法未涉及兩個(gè)顯著性影響因素分析。【擬解決的關(guān)鍵問題】為減小脫粒清選中水稻莖稈EDEM 模擬仿真誤差，筆者以帶有葉鞘包裹的水稻莖稈為研究對象，以莖稈堆積角、莖稈彎曲彈性模量為評價(jià)指標(biāo)，選取EDEM 軟件中HBP（Hertz-Mindlin with bonding particle）模型，結(jié)合Plackett-Burman 試驗(yàn)、Central-Composite 試驗(yàn)和Box-Behnken 試驗(yàn)對水稻莖稈接觸參數(shù)和粘結(jié)參數(shù)進(jìn)行標(biāo)定，對比分析物理試驗(yàn)與仿真試驗(yàn)結(jié)果進(jìn)行驗(yàn)證，以期為水稻植株在脫粒清選等環(huán)節(jié)離散元建模及仿真研究提供參數(shù)依據(jù)。

1 材料與方法

1.1 水稻莖稈本征參數(shù)

選用湖南省常德市桃源縣葉家沖村早稻，品種為‘中早35’。試驗(yàn)在湖南農(nóng)業(yè)大學(xué)物料收獲實(shí)驗(yàn)室完成，隨機(jī)抽10 根水稻，選取整株水稻第一節(jié)（靠近穗頭部分）由葉鞘包裹的莖稈做樣本（圖1），剪去葉片、穗頭與剩余莖稈（水稻植株第一個(gè)結(jié)點(diǎn)之后部分），利用卷尺和數(shù)顯游標(biāo)卡尺測得莖稈本征參數(shù)均值，其長度為201.50 mm、上端（靠近穗頭的一端）外徑為2.69 mm、下端（遠(yuǎn)離穗頭的一端）外徑為3.94 mm、上端壁厚為0.44 mm、下端壁厚1.04 mm，試驗(yàn)莖稈真實(shí)密度為152.3 kg/m3，試驗(yàn)時(shí)莖稈含水率為79.04%。

圖1 試驗(yàn)樣本Fig.1 Test sample

1.2 堆積角試驗(yàn)

依據(jù)文獻(xiàn)[12，15]，采用自制可抽鋼板的水稻莖稈徑向堆積角測定裝置進(jìn)行堆積角試驗(yàn)，鋼板材料為不銹鋼，測定裝置長寬高分別為400，200，200 mm。試驗(yàn)時(shí)所有莖稈長度取200 mm，平均外徑4 mm，平均壁厚0.7 mm，將一定數(shù)量的短莖稈放入堆積角測定裝置一側(cè)，堆積寬度為40 mm，豎直向上迅速抽出鋼板，水稻莖稈整體向另一側(cè)運(yùn)動(dòng)，與水平鋼板形成堆積角，應(yīng)用相機(jī)拍攝記錄莖稈堆積后形成的照片，堆積角測量與處理過程如圖2所示。

應(yīng)用Origin2019軟件對拍攝圖像進(jìn)行圖2過程處理，由擬合直線斜率計(jì)算莖稈堆積角。試驗(yàn)圖像處理時(shí)，為減少拍攝角度帶來的誤差，分別以水平鋼板為X 軸，豎直固定鋼板為Y 軸進(jìn)行曲線擬合，試驗(yàn)重復(fù)5 次，進(jìn)而得到堆積角均值為47.67°，擬合直線決定系數(shù)為0.937 9，接近于1，表明擬合方程可靠程度高，物理試驗(yàn)所測堆積角可用于后續(xù)仿真試驗(yàn)尋優(yōu)及參數(shù)驗(yàn)證。

圖2 堆積角測量與處理Fig.2 Measurement and processing of stacking angle

1.3 三點(diǎn)彎曲試驗(yàn)

筆者采用三點(diǎn)彎曲試驗(yàn)對莖稈彎曲特性進(jìn)行分析測量[17-20]，基于彎曲特性中的彈性模量對莖稈粘結(jié)參數(shù)進(jìn)行標(biāo)定。試驗(yàn)采用CMT6104 電子萬能試驗(yàn)臺，其主要性能參數(shù)如表1 所示。試驗(yàn)時(shí)將莖稈兩端固定在水平支座上，支座間距為30 mm，試驗(yàn)通過萬能力學(xué)試驗(yàn)裝置以0.06 m/s對直徑為6 mm 的鋼制圓柱壓頭從莖稈中心加載（圖3所示），彎曲彈性模量Eb按下式計(jì)算可得[20]：

圖3 三點(diǎn)彎曲試驗(yàn)Fig.3 Three point bending test

表1 CMT6104主要性能參數(shù)Tab.1 Main performance parameters of CMT6104

式中Eb是莖稈彎曲彈性模量，MPa；F是加載力，N；S是兩支撐間的長度，mm；W是莖稈的彎曲撓度，mm；I是莖稈橫截面相對中性軸的慣性矩[20]，mm4。

根據(jù)生物學(xué)特性，水稻莖稈橫截面可近似為圓環(huán)[20]，其慣性矩I按下式計(jì)算：

式中d是莖稈外徑，mm；t是莖稈壁厚，mm。

試驗(yàn)進(jìn)行10次取均值，根據(jù)式 （1）和式（2）可得莖稈彎曲彈性模量為3 740.71 MPa。

2 仿真模型

2.1 莖稈模型創(chuàng)建

筆者采用多球圓環(huán)排列方式構(gòu)建單層水稻莖稈模型[20]，由離散元軟件的顆粒工廠和顆粒體力API顆粒替換構(gòu)建莖稈模型，通過內(nèi)置的HBP 模型沿軸向粘接。由于真實(shí)稈壁很薄，若采用真實(shí)模擬，則計(jì)算量很大，故在模擬中放大了稈壁厚度[20]，依據(jù)文獻(xiàn)[20]，仿真放大了的稈壁厚度對試驗(yàn)結(jié)果無顯著影響，因此，根據(jù)物理試驗(yàn)莖稈基本特征參數(shù)將仿真模型莖稈壁厚設(shè)為1 mm、外徑為4 mm、長度為200 mm。

圖4 水稻莖稈離散元模型Fig.4 Discrete element model of rice stalk

2.2 莖稈-鋼板靜摩擦系數(shù)測定

測定時(shí)將待測莖稈軸向放置（保證莖稈只滑動(dòng)不滾動(dòng)）在帶有三量數(shù)顯傾角儀的鋼板上，三量數(shù)顯傾角儀通過雙面膠固定在鋼板上，鋼板一端與具有60°傾角的固定鋼板平齊，鋼板下放置一滑塊，通過推動(dòng)滑塊以調(diào)節(jié)鋼板的傾斜角度，當(dāng)水稻莖稈與不銹鋼板發(fā)生位移時(shí)停止推動(dòng)，記錄此時(shí)數(shù)顯角度測試儀的讀數(shù)α，按照式（3）計(jì)算摩擦系數(shù)，試驗(yàn)重復(fù)10次取均值，測定過程如圖5所示，測得莖稈-鋼板靜摩擦系數(shù)為0.4。

圖5 靜摩擦系數(shù)測定試驗(yàn)Fig.5 Test for determination of static friction coefficient

式中fs為靜摩擦系數(shù)；α為靜摩擦臨界角，（°）。

2.3 莖稈-鋼板滾動(dòng)摩擦系數(shù)測定

莖稈-不銹鋼板滾動(dòng)摩擦系數(shù)測定方法與靜摩擦系數(shù)測定方法相似，測定時(shí)將莖稈橫向放置在鋼板上，拍攝滑塊前進(jìn)鋼板緩慢抬起的視頻過程，當(dāng)莖稈發(fā)生滾動(dòng)后停止拍攝，用Premiere 軟件導(dǎo)出莖稈開始滾動(dòng)那一幀圖片，讀出此時(shí)數(shù)顯角度測試儀的讀數(shù)β，按照式（7）計(jì)算滾動(dòng)摩擦系數(shù)，試驗(yàn)重復(fù)10次取均值，測定過程如圖6所示，測得莖稈-鋼板滾動(dòng)摩擦系數(shù)為0.3。

圖6 滾動(dòng)摩擦系數(shù)測定試驗(yàn)Fig.6 Rolling friction coefficient measurement test

式中M為滾動(dòng)摩擦力偶矩，N*m；FN2為支撐面的正壓力，N；fr為滾動(dòng)摩擦系數(shù)；G為莖稈的重力，N；β為滾動(dòng)摩擦臨界角，（°）；r為莖稈的半徑，mm。

2.4 仿真參數(shù)選取與模型建立

根據(jù)文獻(xiàn)[21]，試驗(yàn)中不銹鋼板密度為7 860 kg/m3、泊松比為0.3、剪切模量為7.9×1010Pa。在莖稈堆積角仿真試驗(yàn)時(shí)，莖稈生成方式為Dynamic，生成速率設(shè)為5 000 個(gè)/s，生成總數(shù)量為100 個(gè)，時(shí)間步設(shè)為20%，時(shí)間間隔為0.01 s，網(wǎng)格尺寸為3.0R，莖稈生成時(shí)間為0.2 s，待莖稈模型生成后取出隔板進(jìn)行堆積，堆積時(shí)間為2 s。在莖稈三點(diǎn)彎曲仿真試驗(yàn)時(shí)，仿真幾何模型與實(shí)測模型相同，時(shí)間步設(shè)為10%，時(shí)間間隔為0.01 s，網(wǎng)格尺寸為3.0R，莖稈彎曲仿真時(shí)間為2 s，建立仿真模型如圖7所示。

圖7 仿真試驗(yàn)?zāi)Ｐ虵ig.7 Simulation test model

3 試驗(yàn)設(shè)計(jì)與結(jié)果分析

3.1 Plackett-Burman試驗(yàn)及顯著性分析

Plackett-Burman 試驗(yàn)是基于目標(biāo)響應(yīng)與各因素間關(guān)系，比較每個(gè)因素2 水平間的差異來確定因素顯著性[14]。試驗(yàn)以水稻莖稈堆積角θ和彎曲彈性模量Eb為響應(yīng)值，根據(jù)國內(nèi)外文獻(xiàn)對影響參數(shù)設(shè)計(jì)如表2 所示[5-10,13-16]，各參數(shù)取高低2 個(gè)水平，采用1 個(gè)中心點(diǎn)，對仿真接觸參數(shù)和粘結(jié)參數(shù)的顯著性進(jìn)行評估試驗(yàn)。

表2 Plackett-Burman仿真試驗(yàn)參數(shù)列表Table2 Plackett-Burman simulation test parameter list

由Plackett-Burman 仿真試驗(yàn)得到的響應(yīng)值結(jié)果如表3 所示，其中X1～X8為參數(shù)編碼值，利用Design Expert 軟件對響應(yīng)值θ和Eb進(jìn)行方差分析，得到各影響參數(shù)的顯著性如表4 所示。由表4 可知，莖稈-莖稈靜摩擦系數(shù)us與莖稈-莖稈動(dòng)摩擦系數(shù)u（tP＜0.01）對堆積角θ影響極顯著；莖稈-莖稈滾動(dòng)摩擦系數(shù)ut與莖稈切向粘結(jié)剛度kt（b0.01＜P＜0.05）對莖稈彈性模量Eb影響顯著，莖稈粘結(jié)半徑Rb與莖稈法向粘結(jié)剛度kn（bP＜0.01）對莖稈彈性模量Eb影響極顯著。

表3 Plackett-Burman試驗(yàn)結(jié)果Tab.3 Plackett-Burman test scheme and results

表4 Plackett-Burman試驗(yàn)方差分析Tab.4 Analysis of variance of Plackett-Burman test

3.2 Central-Composite試驗(yàn)及回歸模型

為對影響實(shí)測莖稈接觸參數(shù)指標(biāo)θ的us和ut進(jìn)行標(biāo)定，依據(jù)Plackett-Burman試驗(yàn)結(jié)果，以us和ut為影響因素，θ為評價(jià)指標(biāo)，選取3個(gè)中心點(diǎn)進(jìn)行Central-Composite試驗(yàn)，試驗(yàn)因素設(shè)計(jì)與結(jié)果如表5。仿真試驗(yàn)中對θ影響不顯著的其它參數(shù)均取中間值。根據(jù)表5得到關(guān)于堆積角θ的多元回歸擬合方程為：

表5 Central-Composite試驗(yàn)設(shè)計(jì)與結(jié)果Tab.5 Central-composite test design and results

Central-Composite 試驗(yàn)回歸模型方差分析如表6 所示，從表6 可知，該擬合回歸模型P＜0.01，說明堆積角與所得回歸方程有密切關(guān)系；其中各個(gè)參數(shù)的P值均小于0.05、表明了回歸模型的有效性。試驗(yàn)精確度Adeq Precision=8.991 8＞4，說明模型可行且具有良好的精確度；決定系數(shù)R2=0.925 3 趨近于1，表明物理試驗(yàn)與仿真試驗(yàn)擬合度好，試驗(yàn)誤差小。在Desirability 值為1（可靠度最高）的一百組回歸解中尋找與實(shí)測堆積角相近的一組，可得莖稈-莖稈靜摩擦系數(shù)為0.35、莖稈-莖稈滾動(dòng)摩擦系數(shù)0.016。

表6 Central-Composite試驗(yàn)θ方差分析Tab.6 Analysis of variance of Central-composite test θ

3.3 Box-Behnken試驗(yàn)及回歸模型

為對影響實(shí)測莖稈粘結(jié)參數(shù)指標(biāo)Eb的Gb、Rb和knb進(jìn)行標(biāo)定，基于Plackett-Burman 試驗(yàn)因素范圍以及Central-Composite 試驗(yàn)確定的us和ut，依據(jù)Box-Behnken 試驗(yàn)原理，以Rb、knb和ktb為影響因素，Eb為評價(jià)指標(biāo)，選取3 個(gè)中心點(diǎn)進(jìn)行3 因素3 水平Box-Behnken 試驗(yàn)，試驗(yàn)設(shè)計(jì)與結(jié)果如表7。仿真試驗(yàn)中除已標(biāo)定參數(shù)外，對Eb影響不顯著的其它參數(shù)均取中間值。根據(jù)表7，應(yīng)用Design Expert 軟件得到關(guān)于彈性模量Eb的多元回歸擬合方程為：

表7 Box-Behnken試驗(yàn)設(shè)計(jì)與結(jié)果Tab.7 Box-Behnken test design and results

Box-Behnken 試驗(yàn)回歸模型方差分析如表8所示，分析可知，該擬合線性回歸模型P＜0.01，說明彈性模量與所得回歸方程是極顯著的；其中Rb、knb和ktb的P值均小于0.01、說明這3個(gè)參數(shù)對于彈性模量的影響極顯著，表明了回歸模型的有效性。

表8 Box-Behnken試驗(yàn)Eb方差分析Tab.8 Analysis of variance of Box-Behnken test Eb

試驗(yàn)精確度Adeq Precision=28.912 0＞4，說明模型可行且具有良好的精確度；決定系數(shù)R2=0.993 4趨近于1，表明擬合方程可靠度高。在Desirability值為1得一百組回歸解中尋找與實(shí)測彈性模量相近的一組，可得莖稈粘結(jié)半徑Rb為2.331 mm、莖稈法向粘結(jié)剛度knb為4.7×109N/m3、莖稈切向粘結(jié)剛度ktb為3.7×109N/m3。

3.4 標(biāo)定結(jié)果驗(yàn)證

根據(jù)莖稈接觸參數(shù)和粘結(jié)參數(shù)的標(biāo)定結(jié)果，應(yīng)用EDEM 軟件建立最佳參數(shù)離散元仿真模型，以驗(yàn)證標(biāo)定試驗(yàn)結(jié)果的可靠度。試驗(yàn)時(shí)莖稈-莖稈靜摩擦系數(shù)為0.35、莖稈-莖稈滾動(dòng)摩擦系數(shù)為0.016、莖稈粘結(jié)半徑為2.331 mm、莖稈法向粘結(jié)剛度為4.71×109N/m3、莖稈切向粘結(jié)剛度為3.67×109N/m3，莖稈泊松比為0.3，莖稈剪切模量為5.05 GPa，莖稈-莖稈恢復(fù)系數(shù)為0.3，試驗(yàn)重復(fù)5 次，模擬測得莖稈堆積角分別為47.79°、51.59°、43.87°、46.77°、50.86°，平均值為48.27°，與實(shí)際物理堆積角（47.67°）的相對誤差為1.26%，應(yīng)用SPSS 軟件進(jìn)行T檢驗(yàn)，得P1=0.737；模擬測得莖稈彈性模量分別為3 672.21，3 819.74，3 595.82，3 661.68，3 819.74 MPa，平均值為3 713.84 MPa，與實(shí)際物理彈性模量（3 740.71 MPa）得相對誤差為0.72%，應(yīng)用SPSS 軟件進(jìn)行T檢驗(yàn)，得P2=0.584。兩次T檢驗(yàn)所得P值均大于0.05，表明仿真試驗(yàn)與物理試驗(yàn)無顯著性差異。綜上可知，筆者所標(biāo)定的水稻莖稈離散元仿真參數(shù)準(zhǔn)確可靠，可用于水稻離散元仿真建模研究。

4 結(jié)論

筆者基于EDEM 離散元軟件，選用HBP（Hertz-Mindlin with bonding）模型，通過Plackett-Burman 試驗(yàn)、Central-Composite 試驗(yàn)和Box-Behnken 試驗(yàn)對水稻短莖稈樣本的接觸參數(shù)和粘結(jié)參數(shù)進(jìn)行標(biāo)定，主要結(jié)論如下：

（1）根據(jù)Plackett-Burman試驗(yàn)方差分析可知，莖稈-莖稈靜摩擦系數(shù)與莖稈-莖稈滾動(dòng)摩擦系數(shù)對堆積角影響極顯著，莖稈-莖稈滾動(dòng)摩擦系數(shù)與莖稈切向粘結(jié)剛度對莖稈彈性模量影響顯著，莖稈粘結(jié)半徑與莖稈法向粘結(jié)剛度對莖稈彈性模量影響極顯著。

（2）根據(jù)Central-Composite 試驗(yàn)方差分析和回歸方程可得，對莖稈影響顯著的2 個(gè)接觸參數(shù)的一次項(xiàng)、交互項(xiàng)以及二次項(xiàng)對水稻莖稈堆積角影響顯著；根據(jù)Box-Behnken 試驗(yàn)方差分析和回歸方程可得，對莖稈影響顯著的3 個(gè)粘結(jié)參數(shù)的一次項(xiàng)與交互項(xiàng)對水稻莖稈彈性模量影響極顯著，莖稈粘結(jié)半徑的二次項(xiàng)對水稻莖稈彈性模量影響顯著。

（3）以水稻莖稈實(shí)際堆積角與彈性模量為目標(biāo)，對回歸方程進(jìn)行尋優(yōu)求解，得到顯著性接觸參數(shù)的最優(yōu)值為莖稈-莖稈靜摩擦系數(shù)0.35、莖稈-莖稈滾動(dòng)摩擦系數(shù)0.016，顯著性粘結(jié)參數(shù)的最優(yōu)值為莖稈粘結(jié)半徑2.33 mm、莖稈法向粘結(jié)剛度4.71×109N/m3、莖稈切向粘結(jié)剛度3.67×109N/m3。

（4）最優(yōu)接觸參數(shù)得出的仿真試驗(yàn)堆積角與物理試驗(yàn)堆積角的相對誤差為1.26%，基于SPSS對仿真堆積角進(jìn)行T試驗(yàn)，測得P1=0.737＞0.05，證明仿真堆積角與物理堆積角無顯著性差異；最優(yōu)粘結(jié)參數(shù)得出的仿真試驗(yàn)彈性模量與物理試驗(yàn)彈性模量的相對誤差為0.72%，基于SPSS對仿真彈性模量進(jìn)行T試驗(yàn)，測得P2=0.584＞0.05，證明仿真彈性模量與物理彈性無顯著性差異。標(biāo)定結(jié)果驗(yàn)證試驗(yàn)表明，基于HBP 接觸模型所標(biāo)定的接觸參數(shù)與粘結(jié)參數(shù)可用于水稻莖稈離散元仿真，可為后續(xù)水稻整株仿真研究以及聯(lián)合收割機(jī)物料優(yōu)化設(shè)計(jì)提供參數(shù)依據(jù)。

致謝：湖南省高新技術(shù)產(chǎn)業(yè)科技創(chuàng)新引領(lǐng)計(jì)劃項(xiàng)目（2020NK2002）同時(shí)對本研究給予了資助，謹(jǐn)致謝意！