基于EDEM 的水稻殘茬秸稈離散元仿真參數(shù)標定

馬紫濤， 趙智豪， 全偉， 石方剛， 高晨， 吳明亮

（1.湖南農(nóng)業(yè)大學機電工程學院，長沙 410128； 2.湖南農(nóng)業(yè)大學東方科技學院，長沙 410128）

水稻是我國重要的糧食作物之一，其種植面積大、產(chǎn)量高，水稻收獲后秸稈保有量豐富[1-2]。水稻秸稈是一種重要的農(nóng)業(yè)物料資源，可以用作草食動物飼料、生活燃料、農(nóng)作物肥料等，隨著經(jīng)濟的發(fā)展，水稻秸稈越來越多的被用作有機肥的主要原料[3]。水稻秸稈含有豐富的纖維素及氮、磷、鉀、鈣等養(yǎng)分，是有機肥的重要來源，秸稈還田后自然腐爛，不僅可以改善土壤結(jié)構(gòu)、提高土壤蓄水保肥能力，還能促進植物的生長發(fā)育，進而提高作物產(chǎn)量[4-8]。為了使秸稈還田分布均勻、充分利用秸稈有機質(zhì)、減少腐爛周期，研究者提出了先切碎秸稈再進行鋪放的還田方式，并設(shè)計了相應(yīng)的切碎還田機具[9-11]。

切碎還田機具中的切碎裝置對田間秸稈進行切割，需要克服殘茬秸稈自身的連接強度、依靠秸稈之間的摩擦及秸稈與機具內(nèi)壁的摩擦等來實現(xiàn)切割破碎，因此，研究殘茬秸稈在切碎裝置中的運動規(guī)律，有助于優(yōu)化切碎還田裝置各工作參數(shù)。隨著計算機運算能力的提高，離散元法在物料與機具相互作用及切割破碎方面的應(yīng)用得到了全面發(fā)展[12-18]，而建立較為精確的水稻殘茬秸稈模型是分析切碎裝置與秸稈相互作用關(guān)系的前提。馬彥華等[19]采用物理試驗和仿真優(yōu)化設(shè)計相結(jié)合的方法對苜蓿秸稈離散元參數(shù)進行了標定，以物理試驗休止角為目標值，通過正交組合試驗得到了苜蓿秸稈最優(yōu)參數(shù)組合。侯杰等[2]對帶有葉鞘包裹的水稻莖稈進行了離散元仿真參數(shù)標定，以莖稈堆積角和粘結(jié)參數(shù)為評價指標，得到了水稻模型的最佳參數(shù)組合，并通過SPSS軟件對試驗結(jié)果進行了t檢驗。賈洪雷等[11]利用離散元法對稻秸稈之間、稻秸稈與機具之間的相互作用關(guān)系進行了研究，以堆積角為評價指標標定了秸稈與秸稈、秸稈與農(nóng)機具之間的接觸參數(shù)。張濤等[20]利用徑向堆積角相對誤差值為評價指標，應(yīng)用正交試驗方法標定了玉米秸稈與揉碎機錘片、玉米秸稈之間的接觸參數(shù)。已有研究為作物莖稈的離散元建模提供了方法，但未見對水稻收獲后的殘茬秸稈離散元參數(shù)的標定，而殘茬秸稈離散元建模的準確性是秸稈切碎還田裝置仿真優(yōu)化的前提。因此，為了得到準確的離散元仿真參數(shù)，本研究以水稻收獲后的殘茬秸稈為研究對象，以實測的秸稈抗剪力與徑向堆積角為評價指標，分別選取EDEM 軟件中的Hertz-Mindlin with bonding 與Hertz-Mindlin 為接觸模型，利用Box-Behnken 試驗對水稻殘茬秸稈的接觸參數(shù)與粘結(jié)參數(shù)進行標定，通過對比仿真試驗與物理試驗中的秸稈抗剪力、徑向堆積角來進行驗證，以確定水稻殘茬秸稈的離散元仿真參數(shù)，可為秸稈切碎還田裝置在切碎等環(huán)節(jié)的離散元建模及仿真研究提供參考。

1 材料與方法

1.1 試驗材料

以湖南省水稻研究所種植的早稻‘中早35’秸稈為試驗材料。水稻收獲后殘茬秸稈經(jīng)自然風干，測得其平均秸稈長度為（260.50±5.46）mm，隨機抽取10根秸稈，去除秸稈根部1節(jié)，從地面開始選用第3 節(jié)秸稈，利用游標卡尺測得水稻秸稈的平均長度為（80.6±0.53）mm。秸稈的截面形狀近似橢圓形，長軸和短軸取平均值，分別為5.30和4.16 mm，應(yīng)用電子水分測定儀（DHS-10A）測得秸稈含水率為8.80%，測定試驗秸稈的密度為126.4kg·m-3。

1.2 接觸力學參數(shù)的測定

1.2.1靜摩擦系數(shù) 利用水平儀（三量187-211）和鋼板對秸稈與鋼、秸稈與秸稈之間的靜摩擦系數(shù)進行測定。試驗時，將已去除秸稈根部1 節(jié)的秸稈縱向放置在帶有傾角儀的鋼板上，調(diào)節(jié)鋼板的傾斜角度，當緩慢轉(zhuǎn)動到一定角度時水稻秸稈在鋼板上滑動，記錄此時傾角儀上的指示角度，計算得到秸稈與鋼的靜摩擦因數(shù)；選取較圓的秸稈，利用雙面膠將其橫向放置并貼在鋼板上形成秸稈平面，再將需測定的秸稈縱向放置在秸稈平面上，重復(fù)10次取平均值。

1.2.2滾動摩擦系數(shù)測定 滾動摩擦系數(shù)測定儀器和測試方法與靜摩擦系數(shù)相同。將單根秸稈分別橫向放置在鋼板或秸稈平面上，記錄秸稈開始滾動時傾角儀指示的角度，重復(fù)10 次取平均值。

1.2.3碰撞恢復(fù)系數(shù)測定 碰撞恢復(fù)系數(shù)是碰撞前后兩物體沿接觸處法線方向上的分離速度（v0）與接近速度（v1）之比，只與碰撞物體的材料有關(guān)。恢復(fù)系數(shù)的計算公式如下[20]。

式中，H0為試驗秸稈下落時的初始高度，H1為試驗秸稈碰撞后彈起的最大高度與碰撞接觸點的距離。

選擇下落高度為200 mm，試驗過程為：調(diào)整高速攝影儀處于水平位置并與電腦相連，設(shè)置高速拍攝頻率為500 Hz，秸稈從高度尺為200 mm的位置自由下落，與水平面發(fā)生碰撞并彈起一定高度，高速攝影儀記錄相關(guān)數(shù)據(jù)并導(dǎo)入電腦。應(yīng)用Origin 軟件進行數(shù)據(jù)分析，以高度尺與水平面的交點為坐標原點，高度尺方向為Y軸，水平方向為X軸，得到秸稈與鋼板之間的時間-位移曲線圖（圖1），從而計算得到秸稈與鋼之間的恢復(fù)系數(shù)；秸稈與秸稈之間的恢復(fù)系數(shù)測算方法與上述測算方法一致，秸稈下落前在鋼板上平鋪1 層較圓的秸稈，形成秸稈平面，重復(fù)10次取平均值。

圖1 水稻秸稈運動時間-位移關(guān)系Fig. 1 Time-displacement curve of rice straw

1.3 水稻徑向堆積角測定試驗

根據(jù)文獻[21]，采用可抽的鋼板來測定秸稈的徑向堆積角，試驗測定裝置如圖2 所示，鋼板長×寬×高為350 mm×220 mm×150 mm。測定得到秸稈的長度平均值為80.7 mm，莖稈直徑平均值為5.8 mm。將一定數(shù)量的秸稈放入固定擋板與可抽鋼板之間的空間內(nèi)，向上勻速抽出鋼板，水稻秸稈整體向另一側(cè)移動，與水平面形成堆積角，應(yīng)用相機拍攝記錄秸稈堆積后成型的照片（圖2A）。應(yīng)用Matlab 軟件對拍攝圖像進行二值化處理，由擬合直線斜率計算秸稈堆積角（圖2B、C），實測徑向堆積角均值為21.2°。

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

1.4 剪切試驗

采用剪切試驗對水稻秸稈進行力學分析[22-23]，基于水稻秸稈抗剪強度中等的力對秸稈的粘結(jié)參數(shù)進行標定，試驗采用CMT6104 電子萬能試驗臺。試驗時將秸稈兩端放置在夾具上，為防止秸稈發(fā)生滑移，秸稈兩端用螺栓固定，試驗時自制刀片通過萬能力學試驗裝置以100 mm·min-1的速度向下運動并從莖稈中心進行切割。重復(fù)10次取平均值，可得秸稈抗剪力均值為23.8 N。

1.5 秸稈仿真參數(shù)標定

1.5.1秸稈本征參數(shù)設(shè)置 采用離散元仿真軟件EDEM 2021，由離散元軟件的顆粒工廠生成水稻秸稈模型。考慮計算精度與計算量，用8 個直徑為1 mm 的圓球構(gòu)成1 個橢圓環(huán)，設(shè)置橢圓截面長軸為5 mm，短軸為4 mm，用81個這樣的橢圓環(huán)均勻排布形成橢圓柱水稻秸稈，長度為81 mm，仿真計算時顆粒與顆粒之間選用Hertz-Mindlin with bonding 接觸模型，秸稈與秸稈之間選用Hertz-Mindlin 接觸模型，仿真試驗中先進行秸稈剪切試驗，再根據(jù)剪切試驗得到的單個秸稈模型進行徑向堆積角仿真試驗，離散元仿真時水稻秸稈模型見圖3。根據(jù)水稻秸稈的物理本征參數(shù)，結(jié)合相關(guān)文獻[2,20,24]，水稻秸稈離散元仿真本征參數(shù)設(shè)置見表1。

表1 離散元仿真材料本征參數(shù)設(shè)置Table 1 Material of parameter setting of discrete element simulation

圖3 水稻秸稈離散元仿真模型Fig. 3 Discrete element simulation model of rice

1.5.2秸稈接觸參數(shù)獲取與仿真模型建立 為保證真實性，在水稻秸稈剪切試驗時，參照物理試驗在EDEM 中建立固定秸稈的仿真模型，秸稈模型與實際相同，仿真試驗各參數(shù)與物理試驗保持一致，時間步設(shè)為2%，存儲時間間隔為0.01 s，秸稈剪切時間為6 s；在秸稈徑向堆積角仿真試驗時，依據(jù)實測堆積角試驗裝置，按1∶1 建立虛擬仿真模型，秸稈生成方式為Dynamic，生成速率設(shè)為5 000 個·s-1，生成總數(shù)量為80 個，時間步設(shè)為15%，堆積時間為3 s。建立的仿真模型如圖4所示。

圖4 仿真試驗?zāi)Ｐ虵ig. 4 Simulation test model

1.5.3秸稈接觸參數(shù)獲取 離散元秸稈由橢圓球顆粒勻排布形成，與實際秸稈相比粗糙度有一定差異。本研究以實際測定的秸稈接觸參數(shù)為依據(jù)，選擇秸稈與秸稈之間的碰撞恢復(fù)系數(shù)（A）、靜摩擦系數(shù)（B）、動摩擦系數(shù)（C）3 個接觸力學參數(shù)為試驗因素[2]；為了減少Hertz-Mindlin with bonding 接觸模型標定參數(shù)的個數(shù)，已有研究[25-26]表明，該模型條件下粘結(jié)剛度參數(shù)對顆粒的運動規(guī)律影響不顯著，因而可取法向、切向粘結(jié)剛度分別為9×1010、9×106N·m-3。臨界應(yīng)力的取值關(guān)系到顆粒間的強度，進而決定秸稈的抗剪力大小，根據(jù)實測的秸稈外形尺寸、密度和抗剪力大小，結(jié)合已有研究[2,25-26]得到法向、切向臨界應(yīng)力的取值范圍分別為3×1010~7×1010、3×106~7×106Pa。開展表2 所示最陡爬坡試驗，結(jié)果表明3號水平的接觸參數(shù)值所對應(yīng)的秸稈抗剪力相對誤差最小，故各接觸參數(shù)最優(yōu)值的低水平為2號水平值，高水平為4號水平值，即法向臨界應(yīng)力為4×1010~6×1010Pa，切向臨界應(yīng)力為4×106~6×106Pa，粘結(jié)半徑為0.75~1.25 mm。

表2 最陡爬坡試驗結(jié)果Table 2 Design and results of steepest climbing test

選擇粘結(jié)半徑（x1）、法向臨界應(yīng)力（x2）和切向臨界應(yīng)力（x3）3 個接觸參數(shù)作為秸稈剪切仿真試驗的因素，各因素水平范圍值通過最陡爬坡試驗獲得；根據(jù)實測秸稈之間的接觸參數(shù)及已有研究[2,26]，利用Design-Expert 數(shù)據(jù)處理軟件，按照Box-Behnken 試驗設(shè)計方法進行秸稈剪切與徑向堆積角仿真試驗，建立回歸方程與優(yōu)化模型，仿真試驗的因素水平編碼如表3所示。

表3 仿真試驗因素與水平Table 3 Factors and levels of simulation test

2 結(jié)果與分析

2.1 秸稈接觸力學參數(shù)測定結(jié)果

由表4 可知，秸稈與鋼、秸稈與秸稈的靜摩擦系數(shù)分別為0.29、0.23，秸稈與鋼、秸稈與秸稈的滾動摩擦系數(shù)分別為0.17、0.15，秸稈與鋼、秸稈與秸稈之間的碰撞恢復(fù)系數(shù)均值分別為0.22、0.31。

表4 秸稈接觸力學參數(shù)測定結(jié)果Table 4 Results of straw contact mechanical parameters

2.2 秸稈剪切仿真試驗結(jié)果

2.2.1秸稈抗剪力回歸模型與方差分析 根據(jù)試驗因素水平制定二次旋轉(zhuǎn)正交組合試驗方案，進行17組仿真試驗，每組試驗重復(fù)3次，取平均值作為試驗結(jié)果。利用Design-Expert 8.0.6軟件進行試驗方案設(shè)計及結(jié)果分析（表5），構(gòu)建數(shù)學模型，檢驗其顯著性，分析交互作用的影響規(guī)律[27]，對試驗結(jié)果中的數(shù)據(jù)進行二次多元回歸擬合，選用二次項建立抗剪力（F）與粘結(jié)半徑（x1）、法向臨界應(yīng)力（x2）及切向臨界應(yīng)力（x3）之間的回歸模型，并剔除不顯著因素，得到回歸方程如下。

表5 試驗設(shè)計及結(jié)果Table 5 Experimental design and results

由表6 可知，抗剪力模型的P=0.000 2，遠小于0.01，說明回歸方程模型極顯著，具有統(tǒng)計學意義；失擬項不顯著（P=0.215 1>0.05），表明無失擬因素存在，可用該回歸模型代替真實試驗對結(jié)果進行分析。根據(jù)F值可知，各因素對抗剪力的顯著性順序由大到小為法向臨界應(yīng)力、切向臨界應(yīng)力及粘結(jié)半徑；決定系數(shù)R2=0.9667，接近1，變異系數(shù)（CV）與試驗精確度分別為6.7%和11.63，說明該抗剪力擬合回歸模型具有較高的可靠性。

表6 抗剪力方差分析Table 6 Variance analysis of shearing resistance

2.2.2響應(yīng)曲面分析 為確保抗剪力仿真試驗的準確性，以實測秸稈抗剪力為優(yōu)化目標，利用Design-expert8.0 數(shù)據(jù)分析軟件中Optimization 功能進行抗剪力影響因素的優(yōu)化求解，其目標函數(shù)與約束條件如下。

根據(jù)約束條件，求解得到優(yōu)化后的最佳參數(shù)組合為：粘結(jié)半徑1.06 mm、法向臨界應(yīng)力4.77×1010Pa、切向臨界應(yīng)力為4.67×106Pa。由抗剪力回歸模型方差分析（表6）可知，模型線性項x2對抗剪力的影響呈現(xiàn)極顯著性，模型線性項x3對響應(yīng)值的影響也呈現(xiàn)顯著性，模型交互作用項x2x3、二次項對堆積角的影響均呈現(xiàn)顯著性；利用軟件的Numerical 工具進行交互作用分析，設(shè)置粘結(jié)半徑為1.06 mm，考察法向臨界應(yīng)力與切向臨界應(yīng)力的交互作用對秸稈抗剪力的影響，得到響應(yīng)曲面如圖5 所示。由圖5 可知，當粘結(jié)半徑為1.06 mm 時，抗剪力隨著切向臨界應(yīng)力的增大先增大后減小，切向臨界應(yīng)力在4.5×106~5.0×106Pa時抗剪力達到最大值，而隨著法向臨界應(yīng)力的增大抗剪力一直呈增大趨勢。由此可見，當切向臨界應(yīng)力在4.5×106~5×106Pa、法向臨界應(yīng)力在4×1010~6×1010Pa 時存在目標值點，該點即為響應(yīng)曲面圖中的所求點。

圖5 交互作用對抗剪力的影響Fig. 5 Effects of interaction of factors on shearing resistance

2.3 徑向堆積角仿真試驗結(jié)果與分析

2.3.1徑向堆積角回歸模型與方差分析 根據(jù)試驗因素水平制定二次旋轉(zhuǎn)正交組合試驗方案，進行17組仿真試驗，每組試驗重復(fù)3次，取平均值作為試驗結(jié)果。利用Design-Expert 8.0.6 軟件進行試驗方案設(shè)計及結(jié)果分析（表7），構(gòu)建模型，檢驗其顯著性，分析交互作用的影響規(guī)律[27]。對試驗結(jié)果中的數(shù)據(jù)進行二次多元回歸擬合，選用二次項建模建立堆積角與碰撞恢復(fù)系數(shù)（A）、靜摩擦系數(shù)（B）及動摩擦系數(shù)（C）之間的回歸模型，并剔除不顯著因素，得到的回歸方程式如下。

表7 試驗設(shè)計及試驗結(jié)果Table 7 Experimental design and results

由表8可知，堆積角模型的P=0.001 9，遠小于0.01，說明回歸方程模型極顯著，具有統(tǒng)計學意義；失擬項不顯著（P=0.083 1>0.05），表明無失擬因素存在，可用該回歸模型代替真實試驗對結(jié)果進行分析。根據(jù)F值的大小可知，各因素對堆積角的影響順序由大到小為靜摩擦系數(shù)、碰撞恢復(fù)系數(shù)、動摩擦系數(shù)。決定系數(shù)R2=0.938 0，接近1，變異系數(shù)（CV）值與試驗精確度分別為2.68%和11.61，說明該堆積角擬合回歸模型具有較高的可靠性。

表8 堆積角方差分析Table 8 Variance analysis of stacking angle

2.3.2響應(yīng)曲面分析 以實測徑向堆積角為優(yōu)化目標，優(yōu)化求解方法與前述秸稈抗剪力的求解方法一致，得到優(yōu)化后的最佳參數(shù)組合為：碰撞恢復(fù)系數(shù)0.21，靜摩擦系數(shù)0.19，動摩擦系數(shù)0.09。由表8 中堆積角回歸模型方差分析可知，模型線性項B對堆積角的影響呈現(xiàn)極顯著性，模型線性項A、C對響應(yīng)值的影響也呈現(xiàn)顯著性，模型交互作用項BC、二次項A2對堆積角的影響均呈現(xiàn)顯著性；以實測堆積角為目標，利用軟件的Numerical工具進行交互作用分析，設(shè)置恢復(fù)系數(shù)為0.21 mm，考察靜摩擦系數(shù)與動摩擦系數(shù)的交互作用對秸稈徑向堆積角的影響，得到響應(yīng)曲面（圖6）。由圖6 可知，當恢復(fù)系數(shù)為0.21 時，靜摩擦系數(shù)在0.13~0.33 區(qū)域內(nèi)時，徑向堆積角隨著靜摩擦系數(shù)的增大而增大，當動摩擦系數(shù)在0.05~0.25區(qū)域內(nèi)時，徑向堆積角隨著動摩擦系數(shù)的增大呈緩慢增大趨勢。由此可見，當動摩擦系數(shù)、靜摩擦系數(shù)分別在0.13~0.33、0.05~0.25區(qū)域內(nèi)時存在徑向堆積角的目標值點，該點即為響應(yīng)曲面圖中的所求點。

圖6 交互作用對堆積角的影響Fig. 6 Effects of interaction of factors on stacking angle

2.4 優(yōu)化結(jié)果驗證

根據(jù)秸稈剪切試驗與徑向堆積角試驗的優(yōu)化結(jié)果，應(yīng)用EDEM 仿真軟件建立最優(yōu)參數(shù)組合條件下的秸稈模型，驗證標定試驗的可靠性。試驗時設(shè)定秸稈與秸稈之間的碰撞恢復(fù)系數(shù)為0.21，靜摩擦系數(shù)為0.19，動摩擦系數(shù)為0.09，秸稈粘結(jié)半徑為1.06 mm，秸稈法向臨界應(yīng)力為4.77×1010Pa，切向臨界應(yīng)力為4.67×106Pa，試驗重復(fù)5 次。測定的仿真試驗秸稈抗剪力分別為25.1、24.3、25.2、23.5、23.1 N，平均值為24.2 N，與實測抗剪力的相對誤差為1.7%；測定仿真試驗秸稈徑向堆積角分別為22.1°、21.8°、22.3°、21.4°、21.6°，平均值為21.8°，與實測徑向堆積角的相對誤差為2.8%。以上結(jié)果表明，通過標定得到的水稻秸稈參數(shù)準確可靠，可為水稻離散元仿真建模提供參考。

3 討 論

目前，科研人員對作物秸稈進行了相關(guān)建模研究。王洪波等[28]采用物理試驗和仿真優(yōu)化設(shè)計相結(jié)合的方法對玉米秸稈離散元仿真參數(shù)進行標定，為揉碎玉米秸稈的其他仿真試驗提供了參考。陳聲顯[29]通過建立玉米秸稈皮、髓、節(jié)的各向同性力學模型，對秸稈成型機進行了設(shè)計。趙吉坤等[30]以秸稈彎曲試驗為基礎(chǔ)，對水稻秸稈離散元力學模型進行了標定。本研究通過測定水稻秸稈的本征物理特性參數(shù)，基于水稻秸稈剪切與徑向堆積角物理試驗，借助EDEM 仿真軟件對水稻秸稈的接觸參數(shù)和粘結(jié)參數(shù)進行了標定，得到最優(yōu)值為秸稈粘結(jié)半徑1.06 mm，秸稈法向臨界應(yīng)力4.77×1010Pa，切向臨界應(yīng)力4.67×106Pa，秸稈與秸稈之間的碰撞恢復(fù)系數(shù)為0.21，靜摩擦系數(shù)為0.19，動摩擦系數(shù)為0.09。

本研究中，通過對水稻秸稈的離散元仿真建模研究，得到了較為準確的水稻秸稈離散元仿真模型，為秸稈切碎裝置的設(shè)計提供了理論基礎(chǔ)，為割刀與秸稈的相互作用關(guān)系的研究提供了參考。影響離散元仿真標定結(jié)果的因素較多，如秸稈的水分含量、規(guī)格尺寸等，本研究將水稻秸稈視為規(guī)格統(tǒng)一的單元來進行研究，但實際稻田秸稈的外形尺寸有一定的差異，后期應(yīng)進行田間秸稈取樣，得到不同規(guī)格秸稈的分布范圍，建立更接近實際情況的秸稈仿真模型。