基于RSM和NSGA-Ⅱ的苜蓿種子離散元模型參數標定

馬文鵬 尤 泳 王德成 尹世杰 郇曉龍

(中國農業大學工學院， 北京 100083)

0 引言

離散元法(Discrete element method, DEM)在農業裝備研發中應用越來越廣泛,采用離散元法研究散粒體動力學已成為一種發展趨勢[1-2]。在排種器工作過程中，種間相互作用關系復雜，應用離散元法對排種過程進行數值模擬，有助于揭示排種機理，進而優化排種參數和性能[3-6]。紫花苜蓿是一種優質牧草，可以為畜禽提供優質飼料，被譽為“牧草之王”[7-9]。優化苜蓿種子離散元模型和仿真參數，可提高離散元法在苜蓿排種器研究中的準確性。物料物性參數主要包括本征參數和接觸參數，本征參數是物體自身的特性參數(如泊松比、彈性模量和密度等)，顆粒間及顆粒與幾何體間的接觸參數包括碰撞恢復系數、靜摩擦因數和滾動摩擦因數，是兩個物體發生接觸時的物性參數。大多數離散元模型的本征參數與真實參數數值一致，但由于顆粒模型的誤差以及顆粒間復雜的接觸特性，應對部分參數進行重新標定及優化[10-12]。

目前，國內外學者已完成土壤、稻谷、飼料和玉米等離散元模型的建立和仿真參數的標定。于慶旭等[13]基于粘結顆粒模型,以堆積角實測值與仿真值相對誤差為指標，對三七種子離散元模型進行了接觸參數標定。鹿芳媛等[14]利用內部坍塌法、側壁坍塌法和重力平衡法獲取了水稻芽種兩種休止角及滑動摩擦角，通過仿真試驗與實際試驗相結合的方法，對不同含水率水稻芽種離散元的主要接觸參數進行了標定。王云霞等[15]利用兩種接觸材料(有機玻璃板與鋁質圓筒)進行玉米種子堆積角仿真試驗，建立種間靜摩擦因數和種間滾動摩擦因數兩個自變量的二元回歸方程，以實際測量種群堆積角作為已知目標量進行數值求解，求得EDEM中玉米種間摩擦因數。張濤等[16]以堆積角、最大靜摩擦因數以及碰撞恢復系數為指標，標定大豆種子離散元模型與排種盤、攪種輪和有機玻璃間接觸參數。張銳等[17]以實際堆積角為指標，對標準球型和非標準球型沙土顆粒接觸參數進行了標定。彭飛等[18]提出一種基于注入截面法的休止角測定裝置與方法，可通過顆粒堆積體截面的輪廓線直接獲取休止角，并以此進行顆粒飼料離散元模型接觸參數組合的優化。羅帥等[19]提出了通過測定基質含水率預測休止角的方法，并以休止角作為參照對蚯蚓糞基質離散元模型進行參數標定。GHODKI等[20]建立了大豆離散元模型，通過對大豆形狀、靜止角等特征進行比較，獲得了該模型的接觸參數。WANG等[21]將離散元法與物理試驗相結合，以堆積角為指標對不規則形狀玉米顆粒的滾動摩擦因數進行了標定。MOUSAVIRAAD等[22]通過休止角對比試驗，優化了玉米離散元模型的形狀和接觸參數，并利用螺旋輸送機的流量進行了驗證。

目前，在針對離散元模型參數標定的研究中，多以單一的休止角或堆積角為試驗指標，且標定結果僅為一組解。因此，提高標定試驗效率、建立多指標參數優化數學模型，并使用更有效的優化方法對其進行求解，對于獲取更精確的模型參數具有重要的現實意義。本文設計一種可同時測定物料休止角和堆積角的裝置，并提出測試方法。以苜蓿種子為研究對象，進行實際落種與仿真落種試驗，以休止角相對誤差和堆積角相對誤差為試驗指標，通過Plackett-Burman試驗篩選出對指標影響顯著的接觸參數，采用響應面分析法(RSM)建立顯著性參數與各指標之間的二階數學模型，利用非支配排序遺傳算法Ⅱ(NSGA-Ⅱ)進行多目標尋優計算，獲取更具有收斂性和多樣性的Pareto最優解集，并進行排種驗證試驗，以期為苜蓿、黑麥草等小粒牧草種子離散元模型參數標定提供參考。

1 苜蓿種子離散元模型參數標定

休止角和堆積角是表征顆粒物料流動、摩擦等特性的宏觀參數，與接觸材料和物料本身物理特性有關[23]，因此可將其應用于離散元模型的接觸參數標定試驗中。

1.1 實際落種試驗

落種試驗裝置如圖1a所示。該裝置主要由種箱、擋板和圓盤組成，試驗材料為中苜一號紫花苜蓿種子。試驗時，先關閉擋板，在種箱中均勻加入適量種子，約占種箱總容積的2/3，并盡可能保證種子層上表面的水平，將擋板快速打開，苜蓿種子在重力作用下落入圓盤，種箱內兩側對稱形成梯形種堆。休止角θ為種箱上部水平面與種堆斜面所夾銳角；堆積角φ為種箱下部圓盤中水平面與種堆斜面所夾銳角，如圖1b所示。

圖1 實際落種試驗Fig.1 Practical seed-dropping test

為減少人為測量導致的誤差，利用Matlab軟件對采集到的休止角和堆積角圖像依次進行去噪處理、灰度處理以及二值化處理，獲取種堆邊界點，其連線即種堆邊界曲線，最后采用最小二乘法對曲線進行擬合，得到一條直線，直線的斜率即苜蓿種子的實際休止角和堆積角的正切值，如圖2所示,其中左圖為休止角，右圖為堆積角。每組試驗重復5次取平均值(表1)，得到實際試驗中紫花苜蓿種堆休止角和堆積角平均值分別為40.64°和31.28°。

表1 紫花苜蓿實際條件下休止角和堆積角測量結果Tab.1 Measurement results of repose angle and accumulation angle of alfalfa under actual conditions

1.2 落種仿真試驗

根據苜蓿種子的物理特性(表2)，在EDEM軟件中利用球形顆粒組合的方法建立苜蓿種子仿真模型，如圖3所示。顆粒接觸模型選取Hertz-Mindlin無滑移模型[24-25]。將測量裝置的幾何模型和苜蓿種子模型導入EDEM中，設置苜蓿種子數量、試驗過程與真實試驗相同，試驗結束后對堆積角進行圖像采集，同樣利用Matlab處理圖像，獲取仿真休止角β與堆積角γ，如圖4所示。休止角誤差Y1和堆積角誤差Y2的計算公式為

(1)

(2)

表2 中苜一號紫花苜蓿種子物理特性參數Tab.2 Alfalfa seed physics features

圖4 落種仿真試驗Fig.4 Simulated seed dropping test

1.3 Plackett-Burman試驗

應用Design-Expert軟件設計Plackett-Burman試驗，以苜蓿種群休止角誤差和堆積角誤差為響應值，篩選出對響應值影響顯著的仿真參數。仿真試驗中共7個真實參數A～G，設置4個虛擬參數H～K，每個參數按照推薦范圍值均取低、高2個水平，分別以編碼-1和1表示。根據文獻[13-14]以及前期大量預試驗對各參數范圍進行選取，結果如表3所示。仿真試驗中設定1個中心點，共進行12組試驗，每組仿真試驗重復4次，取平均值，并計算休止角誤差和堆積角誤差。

表4為Plackett-Burman試驗的設計方案及仿真結果，采用Design-Expert 軟件對該模擬試驗結果進行方差分析，得出各參數的影響效果如表5所示。

表3 Plackett-Burman 試驗參數Tab.3 Test parameters of Plackett-Burman

由表5可知，種間碰撞恢復系數E、種間靜摩擦因數F和種間滾動摩擦因數G對休止角誤差、堆積角誤差影響顯著，其余參數影響較小。因此在最陡爬坡試驗以及正交試驗中對E、F、G3個影響顯著的接觸參數進行優化。

1.4 最陡爬坡試驗

表6所示為最陡爬坡試驗設計方案及結果，結果表明：隨著種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數的增大，休止角誤差和堆積角誤差呈現先減小后增大的趨勢，取試驗3為中心點，設為中水平，選取試驗2、4水平分別為低、高水平進行三因素二次旋轉正交組合試驗和回歸模型分析；種間碰撞恢復系數、種間靜摩擦因數以及種間滾動摩擦因數的優化范圍分別為0.36～0.54、0.16～0.34、0.050～0.100。

表4 Plackett-Burman試驗方案設計及結果Tab.4 Design and results of Plackett-Burman test scheme

表5 Plackett-Burman 試驗參數顯著性分析Tab.5 Significance analysis of Plackett-Burman test parameters

表6 最陡爬坡試驗設計方案及結果Tab.6 Test design scheme and results of steepest climb

1.5 三因素二次旋轉正交組合仿真試驗

為進一步獲取最優接觸參數組合，以種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數為試驗因素，以休止角誤差和堆積角誤差為試驗指標，進行三因素二次旋轉正交組合仿真試驗。試驗因素編碼見表7，試驗方案與試驗結果見表8(e、f、g為因素編碼值)，每組試驗重復5次取平均值。

表7 仿真試驗因素編碼Tab.7 Simulation test factors and codes

表8 試驗設計方案及響應值結果Tab.8 Test design scheme and response value results

1.6 多目標優化方法

本文選取種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數進行優化，設有兩個響應值，因此屬于多目標優化研究。對于多目標優化問題，通常存在一個最優解集，被稱為Pareto最優解。首先應用RSM法進行正交試驗，建立參數與目標之間的關系，之后利用NSGA-Ⅱ算法進行全局搜索，得到Pareto最優前沿。NSGA-Ⅱ算法是由DEB提出的一種在多目標遺傳算法界非常有效，并基于非支配排序、精英保留策略與多樣性維持機制的優化算法，該算法所特有的精英保留策略和多樣性維持機制可以確保其計算結果的收斂性與多樣性[26-27]。NSGA-Ⅱ算法流程如圖5所示，優化步驟如下：①確定待優化的接觸參數以及試驗指標。②利用RSM進行試驗設計。③從落種仿真試驗中獲取數據。④對試驗結果進行二次方差分析，建立有目標因素的回歸模型。⑤利用NSGA-Ⅱ算法計算Pareto最優前沿，得到最優過程。⑥通過排種試驗驗證最佳工藝參數。

圖5 NSGA-Ⅱ算法流程圖Fig.5 Flow chart of NSGA-Ⅱ algorithm

2 試驗結果與討論

2.1 數學模型的建立與顯著性檢驗

2.1.1休止角誤差回歸模型與顯著性檢驗

通過試驗以及對試驗數據進行多元回歸擬合，得到各因素對休止角誤差影響的回歸模型

Y1=173.771 64-468.323 57E-255.367 7F- 803.419 22G+175.903 72EF-234.130 91EG+ 215.274 73FG+466.093 77E2+354.365 37F2+ 5 472.575 2G2

(3)

回歸方程的顯著性檢驗如表9所示，根據表9可知，該模型的擬合度是極顯著的(P<0.01)。其中E、F、G、EF、E2、F2以及G2的P值均小于0.05，說明以上各項對休止角誤差的影響顯著，相關試驗因素對響應值的影響存在二次關系。對于失擬項P=0.464 3，不顯著，說明不存在其他影響指標的主要因素存在。

2.1.2堆積角誤差回歸模型與顯著性檢驗

通過試驗以及對試驗數據進行多元回歸擬合，得到各因素對堆積角誤差影響的回歸模型

表9 回歸方程方差分析Tab.9 Variance analysis of regression equation

Y2=157.903 25-247.843 26E-515.660 2F- 915.713 49G+356.608 79EF-105.280 34EG+ 42.426 41FG+193.728 02E2+683.234 19F2+ 5 950.715 16G2

(4)

回歸方程的顯著性檢驗如表9所示，根據表9可知，該模型的擬合度極顯著(P<0.01)。其中E、F、G、EF、E2、F2以及G2的P值均小于0.05，說明以上各項對休止角誤差的影響顯著，相關試驗因素對響應值的影響存在二次關系。對于失擬項P=0.064 3，不顯著，說明不存在其他影響指標的主要因素存在。

2.2 各因素對試驗指標的影響

2.2.1各因素交互作用對休止角誤差的影響

對數據進行處理，可得到種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數對休止角誤差的影響，其響應曲面如圖6所示。由圖可知，等高線呈現較大曲率的橢圓形，各因素交互影響顯著。在種間碰撞恢復系數為0.42～0.48、種間靜摩擦因數為0.22～0.28、種間滾動摩擦因數為0.068～0.082時，休止角誤差較小。

2.2.2各因素交互作用對堆積角誤差的影響

各因素對堆積角誤差交互作用的響應曲面如圖7所示。由圖可知，種間靜摩擦因數與種間滾動摩擦因數對堆積角誤差影響較顯著，在種間靜摩擦因數為0.22～0.28、種間滾動摩擦因數為0.068～0.082時，休止角誤差較小。

2.3 NSGA-Ⅱ算法優化

本文以最小休止角誤差和最小堆積角誤差為優化目標，以種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數為優化對象進行研究。在最陡爬坡試驗的基礎上，確定種間碰撞恢復系數為0.36～0.54，種間靜摩擦因數為0.16～0.34，種間滾動摩擦因數為0.05～0.10。因此，優化問題的目標函數和約束函數為

(5)

用NSGA-Ⅱ算法求解優化模型得到的Pareto最優解曲線如圖8所示。對于多目標優化問題，不可能使每個目標同時最優。但可以在目標之間進行協調和權衡，以盡可能地滿足每個目標，這意味著最優邊界上的所有解都可以用于方案優化。在兼顧休止角誤差和堆積角誤差的原則下，選取種間碰撞恢復系數為0.47，種間靜摩擦因數為0.24，種間滾動摩擦因數為0.08，此時休止角誤差為1.687%，堆積角誤差為1.683%。

圖8 基于NSGA-Ⅱ的Pareto最優解曲線Fig.8 Pareto optimal set based on NSGA-Ⅱ

3 驗證試驗

為進一步驗證苜蓿離散元模型和仿真參數的可靠性，采用槽輪式排種器進行試驗驗證，以排種質量流率為試驗指標，在不同排種輪轉速條件下(10～60 r/min)對比質量流率的實測值和仿真值。選取中苜一號苜蓿種子進行臺架試驗，排種輪采用3D打印，臺架如圖9a所示。試驗時，首先在種箱中加入適量紫花苜蓿種子，打開步進電機的電源，調節電機轉速，種子下落至排種器下方的收集盤中，待種子均勻排出時，每隔1 min記錄1次圓盤中種子的質量，計算排種器的質量流率，重復5次取平均值。將排種器的簡化模型、苜蓿種子離散元模型和標定的接觸參數導入離散元軟件EDEM中，在與臺架試驗相同條件下進行仿真試驗，輸出排種器質量流率，仿真試驗如圖9b所示。

圖9 驗證試驗Fig.9 Verification test

圖10為不同轉速條件下質量流率的實測值和仿真值。由圖可知，質量流率的實測值和仿真值與排種輪轉速的變化趨勢一致，兩者平均相對誤差為2.89%，表明該苜蓿離散元模型和接觸參數可用于離散元仿真試驗。

圖10 不同轉速條件下質量流率的實測和仿真結果Fig.10 Measured and simulated values of mass flow rate at different rotational speeds

4 結論

(1)設計了一種可同時獲取物料休止角與堆積角的裝置，以苜蓿種子為研究對象，通過實際落種試驗，并采用Matlab數字圖像處理技術，測定苜蓿種子休止角為40.64°、堆積角為31.28°。

(2)采用Hertz-Mindlin無滑移模型以及球面堆積法，在EDEM軟件中建立苜蓿種子離散元模型，并進行了落種仿真試驗；以休止角相對誤差和堆積角相對誤差為試驗指標，進行Plackett-Burman試驗，篩選出對指標影響顯著的3個因素分別為：種間碰撞恢復系數、種間靜摩擦因數、種間滾動摩擦因數。以這3個因素為試驗因素，以休止角誤差和堆積角誤差為試驗指標，進行了三因素二次旋轉正交組合試驗。通過方差分析可知，3個試驗因素及其二次方項對休止角誤差的影響顯著(P<0.05)。對各因素指標的交互作用進行分析，并獲得了因素與指標之間的回歸數學模型。

(3)以最小休止角誤差和最小堆積角誤差為優化目標，以種間碰撞恢復系數、種間靜摩擦因數和種間滾動摩擦因數為優化對象，采用NSGA-Ⅱ對回歸數學模型進行求解，獲取了Pareto最優解集。在兼顧休止角誤差和堆積角誤差的原則下，選取種間碰撞恢復系數為0.47、種間靜摩擦因數為0.24、種間滾動摩擦因數為0.08，此時休止角誤差為1.687%，堆積角誤差為1.683%。

(4)采用槽輪式排種器進行試驗驗證，搭建排種試驗臺架進行實際排種試驗，并在EDEM軟件中進行仿真試驗。結果表明，在不同排種輪轉速下，苜蓿種子質量流率的實測值和仿真值平均相對誤差為2.89%，該苜蓿種子離散元模型和接觸參數可用于離散元仿真試驗。