飼料油菜薹期收獲莖稈破碎離散元仿真參數標定

廖宜濤 廖慶喜 周 宇 王在騰 蔣亞軍 梁 方

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

0 引言

油菜是我國區域分布最廣、播種面積最大的油料作物，但受農村勞動力轉移、生產過程機械化水平低、種植規模小、進口油料沖擊等多重因素的影響，油菜純油用經濟效益不高[1]。油菜飼料化利用是油菜全價值鏈開發利用的新路徑，主要是在冬春季節刈割抽薹或花期的油菜，全株用于鮮喂或青貯飼喂牲畜，其中薹期油菜鮮嫩多汁，特別適合鮮喂。在冬春季節，飼喂肉牛每日增加5 kg的新鮮油菜，日增重率可提高30%以上[2]；飼喂哺乳母羊，每日添加1 kg的新鮮油菜，能有效促進產后體重恢復，明顯提高羊乳中的脂肪和總固含量[3]；以50%比例油菜，混合稻草、菌渣等為主原料發酵制成的全混合日糧飼喂湖羊，與經青貯玉米為主原料的全日混糧相比，湖羊生長速度和屠宰性狀沒有顯著差異[4]。飼料油菜產量高、營養豐富、成本低，利用水稻、玉米、大豆等作物收獲后的土地空閑期種植，可有效緩解牲畜在冬春季節的青飼料缺乏問題，具有良好的經濟、生態和社會效益，極具發展潛力[5-6]。

飼料油菜生物量大，人工收獲效率低、勞動強度大，機械化收獲是促進飼料油菜快速發展的關鍵。前期研究表明，飼料油菜是含水率高、脆嫩易切斷、易擠壓破損的作物，使用現有飼草收獲機械進行收獲存在過度切碎成糊狀、揉搓擠壓營養流失、輸送通道易堵塞等問題[7]。研制新型飼料油菜收獲機械需要對其切碎、拋送、集料、卸料等多個部件進行結構參數與工作參數的優化設計。近年來，離散元法及其仿真軟件EDEM在農業機械領域得到了廣泛研究與應用[8-13]。應用離散元法研究飼料油菜莖稈與收獲機械的作用機理，可為機具設計與優化提供理論依據，提高研發效率。

飼料油菜的莖稈脆嫩，表面蠟質層、角質層和表皮組織薄，莖稈內外未出現髓腔和木質化，可以用同一屬性顆粒近似構建油菜莖稈物料模型。切碎后的油菜莖稈是散粒體顆粒，在EDEM軟件中可以采用Hertz-Mindlin基本模型；可以通過添加粘結參數Hertz-Mindlin with bonding模型，模擬物料破壞過程[14-17]。離散元仿真模型參數包括材料泊松比、剪切模量和密度等本征參數，碰撞恢復系數、靜摩擦因數和滾動摩擦因數等基本接觸參數，以及模型的法向與切向接觸剛度、法向與切向應力等粘結參數，參數的準確性直接影響仿真結果的可靠性[18-21]。

飼料油菜是一種新型青飼料作物，其物料特性與已有研究物料差異較大，相關研究較少，且缺乏類似物料參考。本文以薹期飼料油菜莖稈為研究對象，以EDEM仿真軟件為平臺，先采用Hertz-Mindlin基本模型，結合飼料油菜莖稈堆積角測試物理試驗，通過因素顯著性篩選和響應曲面試驗，確定碰撞恢復系數、靜摩擦因數和滾動摩擦因數等離散元仿真模型基本接觸參數；再采用Hertz-Mindlin with bonding接觸模型，結合飼料油菜莖稈彎曲破壞試驗，通過響應曲面試驗確定法向與切向接觸剛度、法向與切向應力等粘結參數。以期為薹期飼料油菜機械化收獲過程的離散元仿真分析提供基本參數，并為莖稈類作物切碎過程離散元仿真研究提供模型參數標定方法。

1 材料與方法

1.1 參數測定

1.1.1飼料油菜本征參數

試驗選用華中農業大學現代農業試驗基地種植的“華油雜62”優質雙低油菜，該品種用于飼用，具有產量高、營養指數高的特點[22-23]。種植時采用2BFQ-6型油菜精量聯合直播機播種，播種量5.25 kg/hm2，同步施復合肥675 kg/hm2；平均生長密度46.5株/m2。由于油菜莖稈底部粗壯、頂部細弱，測量統計其莖稈平均直徑為15.18 mm；通過體積計算和質量測定計算薹期莖稈密度為1 049 kg/m3。

選取油菜莖稈植株中間段直徑約15 mm的莖稈，制作長度為20 mm的圓柱顆粒，利用TMS-Pro型質構儀(FTC公司，美國)進行平板單軸壓縮試驗，設置壓縮速度為5 mm/min，加載位移為4 mm，試驗重復10次；通過單軸壓縮試驗前后高度和直徑的變化，計算得到盛花期飼料油菜莖稈彈性模量、剪切模量和泊松比分別為16.61 MPa、5.89 MPa和0.41，變異系數分別為9.65%、9.21%和7.77%。

1.1.2莖稈堆積角

莖稈堆積角通過圓筒提升法試驗[24]。獲得根據飼料油菜鮮喂及青貯的切碎長度要求，將抽薹期飼料油菜的莖稈預制為長度30、40、50、60 mm共4種顆粒樣品，按等質量比例混合均勻；試驗時，將鋼質圓筒(內徑300 mm、高150 mm)固定在萬能材料試驗機上，使底面與試驗臺接觸，向鋼質圓筒內填充飼料油菜顆粒樣品直至填滿，將鋼質圓筒以0.05 m/s的速度向上提升，使樣品形成一個顆粒堆，如圖1所示。試驗重復5次，測量顆粒堆堆積角，統計平均值θ0為33.10°，變異系數為1.27%。

圖1 飼料油菜莖稈顆粒堆積試驗Fig.1 Angle of repose for fodder rape stalk particles

1.1.3莖稈破碎物理參數

為獲得反映粘結模型的法向與切向剛度、法向與切向應力準確參考值，利用TMS-Pro質構儀(FTC公司，美國)以有支撐切割方式慢速對莖稈施加載荷，記錄加載過程彎曲破壞最大載荷作為莖稈破壞力參考值。對薹期飼料油菜莖稈進行預制，獲得長度60 mm、平均直徑15 mm的飼料油菜莖稈進行試驗，固定支撐兩點間距為40 mm，以0.003 m/s速度進行加載，如圖2所示。試驗重復5次，統計最大破壞力平均值為49.94 N，變異系數為4.41%。

圖2 飼料油菜莖稈破壞試驗Fig.2 Bending broken test of fodder rape stalk

1.2 仿真模型

1.2.1莖稈顆粒堆積角仿真模型

在EDEM軟件中，應用Hertz-Mindlin模型進行飼料油菜莖稈堆積仿真試驗。該模型在計算飼料油菜莖稈顆粒與顆粒之間的運動與受力時，實質上是將飼料油菜莖稈顆粒之間的運動與受力過程進行分解，分別為飼料油菜莖稈顆粒間的法向運動、飼料油菜莖稈顆粒間的切向運動和飼料油菜莖稈顆粒間的滾動運動。飼料油菜顆粒近似圓柱體，為匹配實際圓筒提升試驗中飼料油菜莖稈堆積狀態，采用球顆粒組合的方式[25-26]，建立如圖3所示不同長度飼料油菜莖稈模型，為簡化模型，莖稈顆粒直徑取15 mm。結合離散元仿真農業物料和藤莖類物料離散元仿真參數[27-28]，確定本研究中各仿真參數的變化范圍如表1所示。

圖3 飼料油菜莖稈顆粒幾何模型Fig.3 Particle model of fodder rape stalk

表1 莖稈顆粒堆積角仿真模型參數Tab.1 Parameters of angle of repose simulation model for fodder rape stalk particles

注：a表示該項為物理試驗測定，b表示該項為查閱文獻獲得，c表示該項為試驗變量，其范圍為該變量取值的上下限。下同。

在EDEM軟件中，建立與圓筒一致的虛擬圓柱體作為顆粒工廠(內徑300 mm、高150 mm)，顆粒工廠采用動態生成，生成速度為0.4 kg/s，生成時間為4 s，飼料油菜莖稈顆粒總量為1.6 kg。為了保證飼料油菜莖稈顆粒快速穩定，設定顆粒初始下落速度為1 m/s，顆粒生成完成且穩定后，以0.05 m/s的速度向上提升圓筒，仿真顆粒從圓筒底部流出，最終形成一個穩定的顆粒堆，模型如圖4所示。

圖4 莖稈顆粒堆積角仿真模型Fig.4 Simulation model of angle of repose for fodder rape stalk particles

1.2.2莖稈破碎仿真模型

飼料油菜在收獲中涉及到切斷過程，采用Hertz-Mindlin模型建立顆粒模型是一個整體，進行斷裂過程數值計算無效，無法用于飼料油菜切碎過程模擬仿真。在Hertz-Mindlin模型基礎上，加入了顆粒間粘結作用，即采用Hertz-Mindlin with bonding模型，形成具有一定機械強度的粘連顆粒體模型。以該顆粒體模擬生成秸稈，當顆粒間承受的應力達到設置的極限應力時，粘結就被破壞，破壞后的粒子間不再存在粘結力，可用于物料的破碎過程仿真研究。

建模過程如圖5所示。先用小顆粒粘結法生成莖稈模型，替代球形顆粒組合法建立的莖稈模型，具體流程是在EDEM中建立直徑2 mm的小顆粒模型和直徑為15 mm、長度為60 mm、關于原點中心對稱的圓柱體模型，使用小顆粒對圓柱體進行填充[29]；再通過三維軟件Pro/E建立刀具模型，導入EDEM中；最后在EDEM中以刀具為中心建立兩支撐板，支撐板之間間距為40 mm，設定刀具運動速度為0.003 m/s，方向垂直向下，建立莖稈彎曲破壞仿真模型，設置仿真計算固定時間步長為4×10-7s。仿真模型基本接觸參數采用已確定的Hertz-Mindlin模型參數，粘結參數參考藤莖類物料仿真參數[25]，取值范圍如表2所示。

圖5 飼料油菜莖稈彎曲破壞仿真模型Fig.5 Bending broken models of fodder rape stalk

表2 莖稈彎曲破壞仿真模型粘結參數Tab.2 Bonding parameters of bending broken simulation model for fodder rape stalk particle

1.3 參數標定試驗設計

1.3.1基本接觸參數標定

(1)二水平因子試驗設計：應用Design-Expert 10.0.4軟件進行二水平因子試驗設計，篩選顯著影響堆積角的參數。根據表1給定的參數建立莖稈顆粒模型開展堆積角仿真試驗，其中x1～x6分別選擇上、下限兩個數值作為高、低兩個水平；設計16次試驗。

(2)最陡爬坡試驗：針對篩選出的主要參數，進行最陡爬坡試驗，以確定最優值臨近區域。其中非顯著性基本接觸參數取二水平因子試驗中的中間水平，顯著性參數按設定步長逐步增加，建立莖稈顆粒模型開展堆積角仿真分析，記錄仿真試驗得到的堆積角與物理試驗相對誤差變化，根據相對誤差的變化趨勢確定最優值臨近區域。

(3)基本接觸參數響應曲面試驗：基于二水平因子試驗和最陡爬坡試驗結果，根據Box-Behnken Design(BBD)設計原理，取顯著性參數最優值臨近區域的高水平值、中心點和低水平值為響應曲面試驗設計的高、中、低3個水平，開展堆積角仿真試驗。仿真模型中的非顯著性參數取值同最陡爬坡試驗。

1.3.2粘結參數響應曲面試驗

根據Hertz-Mindlin with bonding模型計算原理可知，顆粒間粘結鍵的斷裂與法向接觸剛度、切向接觸剛度、臨界法向應力、臨界切向應力等4個粘結參數相關。為保證參數標定過程中參數范圍的可靠性，避免參數取值超過范圍造成不良影響，根據中心組合設計(Central composite design,CCD)原理，結合表2中參數上下限數值設計響應曲面試驗；以標定的基本接觸參數和設定的粘結參數建立模型，開展飼料油菜莖稈破碎仿真試驗與分析，粘結參數編碼如表3所示。

表3 粘結模型參數編碼Tab.3 Coding of bonding parameters

2 結果與分析

2.1 基本接觸參數標定

2.1.1堆積角影響參數篩選

對影響莖稈顆粒堆積角的基本接觸參數：飼料油菜莖稈間碰撞恢復系數x1、飼料油菜莖稈間靜摩擦因數x2、飼料油菜莖稈間滾動摩擦因數x3、飼料油菜莖稈-鋼碰撞恢復系數x4、飼料油菜莖稈-鋼靜摩擦因數x5、飼料油菜莖稈-鋼滾動摩擦因數x6作為試驗因素，進行二水平部分因子試驗設計，試驗設計與仿真結果如表4所示，方差分析得到各參數的影響效果如表5所示。

表4 二水平因子試驗設計與結果Tab.4 Regular two-level factorial design and results

表5 參數顯著性分析Tab.5 Analysis of parameter significance

由表5可知，根據影響率進行顯著性排序，在飼料油菜莖稈堆積試驗中，x2、x3和x5對堆積角有顯著性影響，其余參數影響極小。因此在后續最陡爬坡試驗和基本接觸參數響應曲面試驗中只考慮這3個影響顯著的參數，其余非顯著性因素取中間水平，即x1取0.5、x4取0.5、x6取0.1。

2.1.2最陡爬坡試驗

由于x2、x3和x5對堆積角的效應均是正值，即堆積角隨參數增加而增加，故設計參數值遞增，觀察并記錄各試驗中堆積角和相對誤差，試驗設計與結果如表6所示。

在3個試驗因素數值逐漸增大的過程中，堆積角隨之增大，與實際物理試驗的堆積角間的相對誤差呈先減后增的趨勢，3號試驗中相對誤差最小，故以3號試驗中的各參數數值作為后期試驗的中心點，2號試驗和4號試驗分別為低水平和高水平進行后續BBD響應面試驗設計。

表6 最陡爬坡試驗設計與結果Tab.6 Experimental results of the steepest ascent test

2.1.3響應曲面試驗結果與分析

根據前面篩選試驗得出的結果，基本接觸參數取x2、x3和x5，利用Design-Expert 10.0.4進行三因素三水平響應曲面試驗設計，中心水平設置5組重復，總共進行17組飼料油菜莖稈堆積角仿真試驗。試驗設計與仿真結果如表7所示。

表7 接觸參數響應曲面試驗設計與結果Tab.7 Response surface test design and results of basic contact parameters

應用Design-Expert 10.0.4軟件對響應曲面試驗結果進行擬合分析，不同數學模型擬合效果對比表明，使用二次全模型方程擬合效果較好，模型決定系數R2為0.974 7，模型顯著(P=0.001 0)且失擬項不顯著，方差分析結果如表8所示，方程為

(1)

由表8可知，二次全模型P小于0.01，且失擬項P大于0.05，證明模型可靠，其中x2、x3和x5以及x3和x5的平方項對堆積角的影響都顯著，但各個因素之間的交互項的P值都比0.05大得多，為了保證響應曲面的生成和堆積角θ的準確性，保留所有因素建立二次多項式回歸方程。

表8 堆積角響應曲面二次全模型方差分析Tab.8 Analysis of variance of angle of repose for response surface quadratic model

通過響應曲面分析所生成堆積角模型，在已知x2、x3和x5的情況下可以確定飼料油菜莖稈顆粒堆積角，但在已知堆積角的情況下，x2、x3和x5的具體數值是堆積角的等值面曲線，因此通過Design-Expert 10.0.4軟件以實際試驗堆積角與仿真試驗堆積角之間相對誤差δθ為響應值開展分析。結果表明二次全模型方程擬合較好，決定系數R2為0.913 7，方程為

(2)

方差分析表明二次多項式回歸模型P為0.000 4，且失擬項P大于0.05，模型可靠，其中除了x2和x5的交互項P遠大于0.05之外，其余因素都具有較好的顯著性。

通過Design-Expert 10.0.4軟件約束求解工具對式(2)求誤差最小極值點，得x2、x3和x5分別為0.43、0.05和0.29。將求解出的參數代入式(1)，可得出堆積角擬合值θ=34.28°，此時計算值與物理實現測定值θ0的相對誤差為3.57%，表明建立的飼料油菜莖稈顆粒堆積角計算回歸模型較好。

為驗證求解的參數的可行性與準確性，利用計算所得的參數建立堆積角仿真模型開展仿真試驗，5次重復試驗測得飼料油菜莖稈平均堆積角為θs=32.35°，與θ0之間的相對誤差δθ=2.27%，模擬結果與實際試驗結果接近，表明確定的飼料油菜顆粒基本接觸參數可行。

2.2 粘結參數標定

利用確定的基本接觸模型參數值和CCD響應面試驗設計確定的法向接觸剛度x7、切向接觸剛度x8、臨界法向應力x9、臨界切向應力x10等4個粘結參數值，開展飼料油菜莖稈彎曲破壞仿真試驗，中心水平設置3組重復，總共進行27組仿真試驗。試驗設計與仿真結果如表9所示。

表9 粘結參數響應曲面試驗設計與結果Tab.9 Response surface test design and results of bonding parameters

對響應曲面試驗結果進行擬合分析，結果表明使用二次全模型方程擬合時，決定系數R2為0.995 4。方差分析(表10)可知，二次全模型P小于0.01，且失擬項P大于0.05；其中法向接觸剛度x7、切向接觸剛度x8及其二者平方項、交互項，臨界法向應力x9和臨界切向應力x10的平方項對彎曲破壞力有顯著影響。

方程影響因子較多，在保證模型顯著、失擬項不顯著的情況下，剔除不顯著項，對二階回歸模型進行優化調整，得到新的回歸方程為

(3)

表10 彎曲破壞力響應曲面二次全模型方差分析Tab.10 Analysis of variance of bending broken force for response surface quadratic model

以破壞力與實測破壞力之間的誤差Δf為響應值進行響應曲面分析，采用二次全模型方程進行擬合，方差分析結果如表11。

表11 彎曲破壞力誤差響應面二次全模型方差分析Tab.11 Analysis of variance of error of bending broken for response surface quadratic model

各因素對破壞力相對誤差Δf的影響顯著性與對破壞力的影響顯著性有較大差異，因素中臨界法向應力x9、臨界切向應力x10以及與這兩項相關的交互項和平方項都不顯著，因此臨界法向應力x9和臨界切向應力x10都取中心水平值，分別為40 MPa和5 MPa。

剔除不顯著項，優化調整后建立破壞力誤差Δf與各顯著性因素之間的二次多項式回歸方程為

(4)

對方程求最小值，得x7=4.60×109N/m，x8=3.55×108N/m。將參數代入式(4)中得Δf=1.61 N，相對誤差為3.22%；將參數代入式(3)中得f=50.95 N，相對誤差為1.98%；表明建立的飼料油菜莖稈顆粒破壞力計算回歸方程可信。以計算所得的參數進行仿真試驗，5次重復平均破壞力為50.64 N，與實際試驗破壞力之間的相對誤差為1.38%，模擬結果與實際試驗結果無明顯差異，標定參數有效。

3 莖稈破碎仿真試驗驗證

為了驗證標定接觸參數和粘結參數的可行性與準確性，確定建立的莖稈破碎離散元模型對于抽薹期飼料油菜莖稈具有普遍適用性，以抽薹期飼料油菜莖稈平均直徑15 mm為中心，分別建立直徑為13、14、16、17 mm的飼料油菜莖稈，模型參數按確定的參數(表12)進行設置，其他設置條件不變，建立離散元模型，與實際試驗對比驗證。

表12 飼料油菜莖稈離散元仿真參數Tab.12 Parameters of fodder rape crop for discrete element simulation

試驗效果如圖6所示，抽薹期不同直徑飼料油菜莖稈破碎狀態無明顯差異，仿真試驗與物理試驗破壞情況一致，表明建立的模型可行。虛擬仿真與物理試驗測定的飼料油菜莖稈破壞力如表13所示。

圖6 莖稈破壞仿真試驗與實際試驗對比Fig.6 Stalk broken simulation experiment and physical experiment

表13 破壞力結果對比Tab.13 Comparison results of broken force

不同直徑飼料油菜莖稈破壞力，仿真試驗和實際試驗結果相對誤差不大于4.21%，且破壞力與莖稈直徑具有良好的線性關系，表明參數標定方法正確，建立的破碎模型準確。

4 結論

(1)試驗測得薹期飼料油菜莖稈平均直徑為15.18 mm，油菜莖稈密度為1 049 kg/m3，莖稈彈性模量、剪切模量和泊松比平均值分別為16.61 MPa、5.89 MPa和0.41；通過物理試驗測得油菜莖稈顆粒的堆積角平均值為33.10°，莖稈彎曲最大破壞力平均值為49.94 N。

(2)飼料油菜莖稈間碰撞恢復系數、飼料油菜莖稈-鋼碰撞恢復系數及飼料油菜莖稈-鋼滾動摩擦因數對堆積角影響較小，其取值確定為0.5、0.5和0.1；飼料油菜莖稈間靜摩擦因數、飼料油菜莖稈間滾動摩擦因數及其平方項、飼料油菜莖稈-鋼靜摩擦因數及其平方項對堆積角的影響均顯著。以實際試驗堆積角與仿真試驗堆積角之間相對誤差為響應值進行分析，建立相對誤差的二次回歸模型，以相對誤差最小為目標，優化求解得出飼料油菜莖稈間靜摩擦因數為0.43、飼料油菜莖稈間滾動摩擦因數為0.05、飼料油菜莖稈-鋼靜摩擦因數為0.29，利用基本接觸參數建立模型，進行堆積角仿真試驗，仿真結果與實測值相對誤差為2.27%。

(3)通過顆粒替代建立飼料油菜莖稈破碎模型，根據CCD原理，建立了粘結參數與破壞力之間的二次回歸模型，由方差分析可知，法向接觸剛度、切向接觸剛度及其二者平方項、交互項，臨界法向應力和臨界切向應力的平方項對破壞力有顯著影響。以破壞力仿真值與實測值之間的誤差為響應值進行分析，建立二次回歸模型，優化求解得到法向接觸剛度、切向接觸剛度、臨界法向應力和臨界切向應力分別為4.60×109N/m、3.55×108N/m、40 MPa和5 MPa。利用粘結參數進行仿真試驗，與實際破壞力之間的相對誤差為1.38%。

(4)以標定的參數建立不同直徑飼料油菜莖稈彎曲破壞仿真模型，仿真試驗與實際試驗破壞情況一致，兩者相對誤差不大于4.21%，且破壞力與莖稈直徑具有良好的線性關系，說明標定的參數準確、可靠，建立的破碎模型正確、可行。

本研究以物理試驗為參考，通過堆積角試驗和對莖稈緩慢施加載荷的彎曲破壞試驗，確定了飼料油菜薹期收獲莖稈破碎離散元模型參數，為建立莖稈類作物切碎過程離散元仿真模型提供了一種方法。受油菜莖稈生物體物料的復雜性影響，莖稈彎曲破壞試驗的取樣僅以抽薹期莖稈中間段、直徑約為15 mm的物料為代表，標定參數具有一定的適用范圍。后續將針對不同生長期和不同部位莖稈進行系統性研究，并結合飼料油菜收獲仿真分析，對模型參數進行修正，對飼料油菜收獲切碎、輸送關鍵裝置進行優化改進。