基于堆積試驗的玉米包衣種子離散元參數標定

李飛翔， 王鵬， 王云飛， 葛越鋒， 唐凱懌， 李得志

（洛陽智能農業裝備研究院有限公司，河南 洛陽 471000）

氣力式排種器因具有排種精度高、對種子損傷少等優點而被廣泛應用。在實際排種作業中，種子處于氣固兩相環境中，運動情況較為復雜，離散元（discrete element method，DEM）與計算流體動力學（computational fluid dynamics，CFD）耦合法能夠有效分析排種器中流場的變化以及種子運動的受力情況，為排種器的結構設計及參數優化提供依據[1-3]。為提高玉米的產量及質量，常常對玉米種子進行包衣處理。種子包衣是將殺蟲劑、微生物肥料、著色劑等非種子材料包裹在種子表面，以提高種子抗逆性、抗病性，對農業機械化、現代化及農業可持續發展有著重要作用。研究人員基于離散元法對包衣種子已做了大量研究：閆聰杰[4]設計一種大豆包衣裝置，并基于離散元法對其包衣過程進行分析；曲芳[5]針對現有大豆種子包衣攪拌裝置包衣不均、破碎率高的問題，基于離散元法設計一種攪拌葉片并優化設計；胡建平等[6]以番茄磁粉包衣種子為對象，基于離散元法研究了種箱振動頻率、振幅對種群運動規律及種箱供種性能的影響。

DEM-CFD 常采用 EDEM 與 Fluent 軟件的耦合，要求EDEM 中顆粒體積小于Fluent 中網格體積。因此，EDEM 中常采用黏結顆粒模型（bonded particle model，BPM），用較小的顆粒黏結成物料顆粒模型[7]。BPM 需設置單位面積法向與切向剛度、極限法向與切向應力、黏接半徑等參數[8]，現階段對于玉米種子的仿真參數設置研究較少。基于離散元的Hertz-Mindlin 模型以堆積角為響應，王云霞等[9]以玉米顆粒種間靜摩擦與滾動摩擦系數為因素，設計試驗對其進行標定；王美美等[10]對玉米籽粒離散元仿真參數進行標定，為玉米籽粒離散元仿真提供參考；李秀清[11]基于離散元法對玉米粉碎過程進行分析，采用測定的楔形玉米籽粒建模，借鑒玉米籽粒壓縮試驗數據，將極限法向應力、極限切向應力作為固定值，選取法向剛度、切向剛度、黏接半徑3個因素進行標定，得到玉米籽粒的仿真接觸參數。

為得到玉米包衣種子BPM 離散元仿真所需的精確的種間接觸參數，本研究將利用激光掃描儀得到的玉米包衣種子輪廓作為顆粒模板，采用BPM 以小顆粒黏結方式生成玉米包衣種子元顆粒模型，以堆積角作為響應值，基于Isight 軟件的近似模型與DOE 聯合模塊，通過因素顯著性篩選、最優逼近以及響應面優化對玉米包衣種子的仿真接觸參數進行標定，以期為玉米包衣種子氣力式排種器的兩相耦合仿真提供參考。

1 材料與方法

1.1 試驗材料獲取

供試玉米品種為鄭單958，純度＞96%，凈度＞90%，發芽率＞90%，含水率＞13.5%，由河南金博士種業股份有限公司提供。由于種子形狀多樣，將鄭單958 玉米種子分為大扁形、長扁形、類球形3類（圖1），比例為2∶3∶5。

圖1 玉米包衣種子形狀Fig.1 Shape of corn coated seed

1.2 接觸模型

為方便與Fluent 軟件耦合，EDEM 顆粒間接觸選用BPM。該模型接觸情況下，黏結鍵可以抵抗切向和法向運動，直至達到最大的法向和切向極限應力為止。黏結作用產生后，元顆粒所受的力和力矩設置為0，并根據每個時間步長進行增量調整[12-13]。

式中，δFn為法向力增量（N）；δFτ為切向力增量（N）；δMn為法向力矩增量（N·m）；δMτ為切向力矩增量（N·m）；Y為基礎區域面積（m2）；Rb為黏結半徑（m）；J為截面極慣性矩；Sn為法向剛度（N·m-3）；Sτ為切向剛度（N·m-3）；vn為法向速度（m·s-1）；vτ為切向速度（m·s-1）；wn為法向角速度（r·s-1）；wτ為切向角速度（r·s-1）；δt為時間步長（s）。

當法向和切向剪應力超過預定值σmax和τmax時，黏結就會斷裂。

黏結應力與剛度值的不同設置代表了不同的材料屬性，高的剛度值將產生高的結合力和應力。小顆粒在黏結之前，根據Hertz-Mindlin 模型相互作用，當達到鍵形成時間時，被定義的小顆粒將會黏結在一起形成元顆粒。

1.3 堆積角測定方法

依據已有對堆積角的研究[14-18]，以漏斗法對玉米包衣種子的堆積角進行測定，如圖2所示，測定5次取其平均值，測得鄭單958玉米包衣種子的堆積角為29.56o±0.27o。由于標定參數為顆粒間接觸，為消除顆粒與部件的影響，采用亞克力材質加工落料底盤，底盤直徑（L）為150 mm，底盤凹面深度為30 mm。待底盤落滿種子后，繼續落種；待玉米包衣種子溢出底盤后，觀察落料中顆粒堆積高度的變化；待堆積高度無顯著性變化時，用鋼尺測定包衣種子的堆積高度（h），由式（9）計算玉米種子堆積角。

圖2 堆積試驗Fig.2 Stacking test

1.4 仿真試驗

玉米種子為非規則形狀顆粒，為精確建立種子輪廓模型以提高標定精準性，選擇不同種類形狀的種子，采用海克斯康Tango-S手持式激光掃描儀對種子的外輪廓進行掃描。得到種子的點云數據后，通過CATIA 軟件對種子點云數據進行處理，最終得到玉米包衣種子的輪廓模型。

仿真模型采用與實際試驗模型相同尺寸，按照篩選的大扁形、長扁形、類球形種子2∶3∶5的比例生成顆粒。顆粒堆積仿真試驗如圖3 所示，為節省仿真時間，先以圓柱底盤為虛擬工廠生成顆粒，待底盤表面落滿顆粒后，導出仿真時間為0時的文件；再以漏斗為虛擬工廠，進行堆積角仿真試驗，待底盤溢出一定數量顆粒后，停止生成顆粒。采用后處理自帶的量角器，測量堆積角。

圖3 玉米種子顆粒堆積的模擬仿真Fig.3 Simulation test model for repose angle

參照玉米種子BPM 參數的設置[19-22]，本研究各參數取值范圍如表1 所示。仿真本征參數設定為玉米種子泊松比0.4，密度1 197 kg·m-3，剪切模量1.37×108Pa[23-26]。

表1 離散元仿真參數Table 1 Parameter in DEM simulation

1.5 標定方法

以玉米包衣種子實際堆積角為響應值，基于Isight 軟件的近似模型和DOE 聯合模塊，通過因素顯著性篩選、最優逼近以及響應面優化對玉米包衣種子的仿真接觸參數進行標定。本研究的Plackett-Burman 設計以玉米包衣種子堆積角為響應值對BPM 參數的顯著性進行篩選。低水平設定為最初原始水平，高水平設為低水平的2倍，為方便對試驗誤差進行分析，同時設定虛擬參數，試驗參數如表2所示。利用Isight軟件的近似模型與DOE模塊聯合對Plackett-Burman結果進行方差分析，運行前，添加ANOVA Table、主效應圖、Pareto圖對各參數的顯著性進行分析。Plackett-Burman試驗后，以對玉米包衣種子堆積角影響顯著的參數進行最陡爬坡，以便能快速進入到最優值的附近區域。最陡爬坡試驗根據Plackett-Burman試驗所得的回歸系數來確定爬坡步長。

2 結果與分析

2.1 輪廓模型掃描分析

圖4 為玉米包衣種子掃描的輪廓模型，可以看出，通過點云處理后得到的不同種類種子輪廓與實際種子輪廓一致，能夠用于仿真分析。

2.2 種子模型填充分析

通過激光掃描得到種子輪廓后，導入EDEM軟件中作為顆粒模板。結合玉米氣力式排種器研究進展[27-28]，此處以0.5 mm 的小球顆粒對玉米種子顆粒模板進行填充，如圖5 所示。對顆粒模板進行上半部透明化顯示，填充的顆粒與模板貼合程度良好，通過小顆粒間形成黏結鍵建立玉米種子元顆粒黏結模型，建立的顆粒模型與玉米包衣種子外形相似。

圖5 玉米包衣種子仿真模型Fig.5 Simulation model of corn coated seeds

2.3 參數標定結果分析

2.3.1 Plackett-Burman 試驗結果分析 Plackett-Burman 設計及結果如表3 所示，堆積角范圍在23.03°～30.31°之間，各仿真值沒有特別跳躍性情況。

表3 Plackett-Burman 試驗設計及結果Table 3 Design and results of Plackett-Burman test

2.3.2 Plackett-Burman 試驗顯著性分析 主效應、Pareto結果如圖6所示，可清晰觀察各個因素對響應的影響趨勢；ANOVA結果如表4所示，可知對玉米包衣種子堆積角影響極其顯著的參數為：玉米種子-玉米種子靜摩擦系數、法向剛度、切向剛度。其余參數結合已有研究進展[29-30]取值（玉米種子-玉米種子恢復系數0.18、玉米種子-玉米種子滾動摩擦系數 0.05、法向應力 1.1×107N·m-3、切向應力4.1×106N·m-3、黏結半徑0.7 mm）進行最陡爬坡以及響應面試驗設計。

表4 Plackett-Burman 試驗參數顯著性分析Table 4 Analysis of significance of parameters in Plackett-Burman test

圖6 玉米種子對顆粒堆積響應分析Fig.6 Response analysis of corn seed to grain accumulation

2.3.3 最陡爬坡試驗結果分析 最陡爬坡試驗設計及結果如表5所示，2號水平仿真堆積角與試驗值的誤差最小。由此選取2 號水平為中心點，1號、3號水平為低、高水平進行后續優化設計。

表5 最陡爬坡試驗及結果Table 5 Design and results of steepest ascent test

2.3.4 Box-Behnken 試驗結果分析 根據最陡爬坡試驗結果，選取顯著性參數進行3個水平試驗設計。選取3個中心點對誤差進行評估，結果如表 6 所示，基于 Isight 軟件 RSM 優化模塊，建立3個參數與堆積角間的二階回歸方程。

表6 Box-Behnken 試驗及結果Table 6 Design and results of Box-Behnken test

Box-Behnken 試驗方差分析結果如表7 所示，玉米包衣種間靜摩擦系數(H)、法向剛度(J)、切向剛度（K）、法向剛度-切向剛度（JK）、靜摩擦系數-法向剛度（HJ）、靜摩擦系數二次項（H2）以及切向剛度二次項（K2）的P值＜0.05，表明這些參數對堆積角的影響顯著。失擬項P=0.093＞0.05，變異系數CV=0.801%，表明模型良好，試驗有較高的可靠性。決定系數R2=0.997 4，校正決定系數RAdj2=0.992 8，預測決定系數RPre2=0.958 9，表明模型能真實地反映實際情況。試驗精密度Adep precision=46.973，表明模型良好的精確度。

表7 Box-Behnken 試驗方差分析Table 7 ANOVA of quadraticl of Box-Behnken test

2.3.5 回歸模型交互效應分析 根據優化回歸模型方差分析結果可知，種間靜摩擦系數-法向剛度（HJ）以及法向剛度-切向剛度（JK）交互項對玉米包衣種子堆積角影響顯著（P＜0.05）。在玉米種間靜摩擦系數（H）為0.2 以及切向剛度（K）為0.5 情況下，應用Isight 軟件對靜摩擦系數-法向剛度（HJ）以及法向剛度-切向剛度（JK）交互效應的響應曲面進行繪制，如圖7所示，可以直觀反映交互項對堆積角的影響。由HJ曲面可知，相對玉米包衣種間靜摩擦系數（H），法向剛度（J）的效應面曲線比較陡，表明其對堆積角影響較為顯著。由JK曲面可知，相對于法向剛度（J），切向剛度（K）的效應面曲線比較陡，表明其對堆積角影響較為顯著。

圖7 HJ與JK的交互效應Fig.7 Interaction effect of HJ and JK

2.3.6 試驗優化結果分析 應用Isight軟件RSM優化模塊，以玉米包衣種子實際測定的堆積角29.56o為目標，對回歸方程尋優求解得到玉米包衣種間靜摩擦系數為0.269，法向剛度為2.54×108N·m-3，切向剛度為5.93×107N·m-3。用所得最佳參數組合進行堆積仿真試驗，結果如圖8 所示。仿真試驗所得堆積角為29.85o，與實際試驗值29.56o誤差為0.98%，表明最佳參數組合可行。

圖8 仿真與實際試驗對比Fig.8 Comparison between simulation and actual test

3 討論

現有堆積試驗多采取物料直接堆積于材料表面，實際仿真標定中，物料與部件間的參數設置對堆積結果影響較大，為減小試驗誤差，堆積試驗采用凹槽底盤，使得物料落滿底盤后，在物料表面形成堆積。本研究選用鄭單958 玉米包衣種子，在對種子進行分類篩選后，通過激光掃描儀對輪廓較好的種子進行掃描得到點云數據，采用逆向工程技術，通過CATIA 軟件對點云數據處理，最終得到玉米種子仿真模型。

根據主效應和Pareto 圖分析，篩選出對堆積角影響顯著的參數：玉米種子-玉米種子靜摩擦系數、法向剛度、切向剛度。基于RSM 模塊，建立3個顯著性參數與堆積角間的二次回歸模型，由方差分析結果可知，交互項靜摩擦系數-法向剛度（HJ）、法向剛度-切向剛度（JK）以及3個顯著性參數的二次項（H2）對玉米包衣種子堆積角影響極顯著。

以玉米包衣種子實際測定堆積角為目標，對回歸方程進行尋優求解，得到參數的最佳組合為：種間靜摩擦系數0.269，法向剛度2.54×108N·m-3，切向剛度5.93×107N·m-3。用所得最佳參數組合進行仿真試驗，與實際值誤差為0.98%，仿真與實際試驗值差異不顯著，表明標定所得參數可用于玉米包衣種子離散元BPM 仿真，為玉米氣力式排種器的耦合模擬提供參考。