基于X射線斷層掃描的小麥籽粒模型構(gòu)建及離散元參數(shù)標(biāo)定

朱 康， 惠延波， 周 穎， 王宏曉， 張?jiān)讫垼?李 輝， 陳艷雷， 鄭金鋒， 白路遙

(河南工業(yè)大學(xué)先進(jìn)制造研究所，鄭州 450001)

離散元法又稱離散單元法，是一種用于模擬并分析散體介質(zhì)系統(tǒng)動力學(xué)行為的數(shù)值方法，廣泛應(yīng)用于地質(zhì)工程、土木工程、農(nóng)業(yè)工程[1]。在農(nóng)業(yè)工程中，農(nóng)業(yè)裝備會涉及到大量的谷物顆粒的運(yùn)動，運(yùn)用離散元法研究兩者之間的相互作用關(guān)系，可以優(yōu)化農(nóng)業(yè)裝備的機(jī)械結(jié)構(gòu)、提升裝備工作效率，為現(xiàn)代農(nóng)業(yè)裝備的數(shù)值化設(shè)計(jì)提供理論基礎(chǔ)[4]。

在進(jìn)行離散元仿真時，模型的精度和仿真參數(shù)的設(shè)置會對仿真結(jié)果產(chǎn)生較大的影響[11]。劉凡一等[9]使用游標(biāo)卡尺測量小麥籽粒的平均長軸、短軸，將小麥近似地看作規(guī)則的橢球型，在EDEM中采用5球組合的方式構(gòu)建小麥離散元模型對小麥離散元參數(shù)進(jìn)行了標(biāo)定，但在實(shí)際中，小麥籽粒形狀是不規(guī)則的，通過數(shù)量較少的球形顆粒的堆疊，無法準(zhǔn)確地反映小麥籽粒的真實(shí)輪廓。趙武云[10]對玉米種子進(jìn)行切片，每片厚度約為1 mm，采集玉米種子剖面圖像，根據(jù)剖面圖像提取種子外輪廓樣條曲線，在SolidWorks中進(jìn)行輪廓曲線放樣繪制玉米種子三維模型，但切片數(shù)量過少，無法完整地表示玉米的輪廓特征。張榮芳等[11]借助三維掃描和逆向工程的方法獲取水稻種子外形輪廓，構(gòu)建不同半徑的填充球離散元模型對種間靜摩擦因數(shù)和動摩擦因數(shù)進(jìn)行了標(biāo)定，但三維激光掃描儀針對尺寸較小的物體掃描的精度較低，無法完整地獲取物體三維輪廓。

因此，本實(shí)驗(yàn)以小麥為研究對象，基于X射線斷層掃描技術(shù)獲取1 024張不同灰度值的二維切片，采用圖像濾波、圖像分割等算法構(gòu)建出小麥籽粒的三維外輪廓模型；在EDEM中導(dǎo)入外輪廓模型，分別填充4種不同數(shù)量的球(17、49、80、126球)來擬合小麥籽粒的真實(shí)形狀，構(gòu)建了更加精準(zhǔn)的小麥籽粒離散元模型；通過休止角實(shí)際實(shí)驗(yàn)與EDEM仿真實(shí)驗(yàn)，運(yùn)用Plackett-Burman實(shí)驗(yàn)、最陡爬坡實(shí)驗(yàn)、響應(yīng)面優(yōu)化實(shí)驗(yàn)標(biāo)定了小麥離散元仿真參數(shù)；最后對標(biāo)定的參數(shù)進(jìn)行仿真分析并與實(shí)際實(shí)驗(yàn)進(jìn)行對比，驗(yàn)證了本方法構(gòu)建的離散元模型標(biāo)定參數(shù)的準(zhǔn)確性，為其他不規(guī)則谷物模型的構(gòu)建提供了參考。

1 小麥籽粒模型構(gòu)建

1.1 實(shí)驗(yàn)材料

實(shí)驗(yàn)材料選取豐德存麥5號，含水率12.5%。

1.2 圖像采集

采用YXLON Y.Cheetah型微米X射線斷層掃描系統(tǒng)(Micro-CT)采集小麥籽粒三維灰度圖像，如圖1所示。隨機(jī)選取外形尺寸完好的小麥籽粒，通過X射線斷層掃描系統(tǒng)對其進(jìn)行掃描，單顆籽粒1次掃描約20 min，可以無損地獲取小麥籽粒的1 024張二維灰度圖像。

圖1 二維灰度圖像

1.3 圖像重構(gòu)

從二維灰度圖像可以看出，小麥籽粒的不同組織對X射線的吸收不同，呈現(xiàn)的灰度值也就不同，因此可以有效地區(qū)分各組織的輪廓[12]。通過二維灰度圖像的堆疊，可完整地構(gòu)建出小麥籽粒輪廓。采用非局部均值濾波算法，去除圖像噪聲，提高圖像質(zhì)量；采用交互式閥值分割算法分離出小麥籽粒的外輪廓組織，如圖2所示；最后基于體繪制及面繪制方法重構(gòu)出小麥籽粒三維外輪廓模型。構(gòu)建好的三維模型存在表面粗糙、孔洞等問題，將三維模型轉(zhuǎn)化為點(diǎn)云格式導(dǎo)入Geomagic Design中進(jìn)行孔洞的填補(bǔ)、表面光滑處理等操作，最終生成實(shí)體模型[13]，如圖3所示。

圖2 圖像處理

圖3 小麥籽粒模型構(gòu)建

1.4 離散元模型構(gòu)建

將X射線斷層掃描方法構(gòu)建的小麥籽粒外輪廓模型導(dǎo)入EDEM中，以此為基礎(chǔ)擬合小麥籽粒的真實(shí)形狀，構(gòu)建出17、49、80、126球組合的小麥籽粒仿真模型，如圖4所示。

圖4 小麥籽粒離散元仿真模型

2 離散元仿真參數(shù)標(biāo)定

2.1 休止角實(shí)際實(shí)驗(yàn)

休止角是物料由高處滑落至平面上堆積形成的錐體的斜面與水平面之間的夾角，是表征顆粒流動、摩擦特性等特性的宏觀參數(shù)[14]。運(yùn)用離散元法研究的谷物的休止角，可以掌握谷物堆積形成糧堆的機(jī)制，了解顆粒與顆粒、顆粒與裝備的之間的力學(xué)特性。

本次離散元參數(shù)標(biāo)定采用休止角實(shí)驗(yàn)與離散元仿真實(shí)驗(yàn)結(jié)合的方法，實(shí)驗(yàn)方法參照GB/T 16913—2008，實(shí)驗(yàn)裝置如圖5所示，由無底圓筒(長度120 mm，內(nèi)徑40 mm)、有機(jī)玻璃板(150 mm×150 mm)、實(shí)驗(yàn)臺組成。無底圓筒與底面有機(jī)玻璃板重合，將小麥填滿無底圓筒后，以0.02 m/s的速度向上勻速提升圓筒，小麥自然落下并在底面堆積，錐面與底面之間的夾角即為休止角，每組實(shí)驗(yàn)重復(fù)5次。采集籽粒堆積的8個方向的圖像，使用Matable讀取采集的圖像，對圖像進(jìn)行二值化、邊緣提取、邊緣直線擬合來測量小麥堆積的休止角[15]，如圖6所示。對采集的數(shù)據(jù)去除極值，取其平均值，最終測得實(shí)際休止角為20.326 5°。

注：1.實(shí)驗(yàn)臺 2.小麥 3.無底圓筒 4.有機(jī)玻璃板。圖5 小麥休止角實(shí)驗(yàn)

圖6 Matable圖像處理

2.2 休止角仿真實(shí)驗(yàn)

2.2.1 仿真參數(shù)

本次實(shí)驗(yàn)所需的仿真參數(shù)如表1所示[16]。

表1 小麥離散元仿真參數(shù)

2.2.2 仿真實(shí)驗(yàn)

為縮短仿真時間，休止角仿真實(shí)驗(yàn)選用17球小麥模型，在EDEM中構(gòu)建長120 mm、內(nèi)徑40 mm、上下開口的圓筒。在圓筒頂端以5 000顆/s的速度隨機(jī)生成小麥顆粒直至填滿圓筒，等待顆粒穩(wěn)定，此過程時間為1 s；隨后以0.02 m/s的速度向上提升圓筒，直至顆粒全部落下并形成穩(wěn)定的堆積體，此過程過程仿真時間為3 s，如圖7所示。仿真模型選用Hertz-Mindlin (no slip)接觸模型，時間步長為20%，仿真總時長為5 s。

圖7 小麥籽粒休止角離散元仿真

2.3 參數(shù)標(biāo)定

2.3.1 Plackett-Burman實(shí)驗(yàn)

在Plackett-Burman實(shí)驗(yàn)設(shè)計(jì)的基礎(chǔ)上，以小麥籽粒的休止角為響應(yīng)值，選取影響小麥籽粒休止角的八個參數(shù)進(jìn)行評價，篩選出影響小麥休止角的顯著因素[19]。每個因素取低水平(-1)和高水平(+1)2個水平，共12組實(shí)驗(yàn)，如表2所示。

表2 Plackett-Burman實(shí)驗(yàn)設(shè)計(jì)的因素水平及編碼

Plackett-Burman實(shí)驗(yàn)設(shè)計(jì)及結(jié)果如表3所示。運(yùn)用Design-Expert軟件進(jìn)行求解得到的結(jié)果如表4所示。根據(jù)各因素的P值大小判斷出影響小麥休止角的關(guān)鍵因素依次為：E(小麥-小麥靜摩擦系數(shù))>G(小麥-小麥滾動摩擦系數(shù))>F(小麥-有機(jī)板靜摩擦系數(shù))。其中小麥-小麥靜摩擦系數(shù)、小麥-小麥滾動摩擦系數(shù)及小麥-有機(jī)板靜摩擦系數(shù)均達(dá)到極其顯著水平。在后續(xù)的最陡爬坡實(shí)驗(yàn)中只考慮這3個參數(shù)的影響，結(jié)合其他仿真參數(shù)進(jìn)行最陡爬坡實(shí)驗(yàn)和響應(yīng)面優(yōu)化實(shí)驗(yàn)。

表3 Plackett-Burman實(shí)驗(yàn)設(shè)計(jì)及結(jié)果

表4 Plackett-Burman實(shí)驗(yàn)設(shè)計(jì)各因素效應(yīng)評價

2.3.2 最陡爬坡實(shí)驗(yàn)

基于Plackett-Burman實(shí)驗(yàn)篩選出較顯著的參數(shù)，設(shè)計(jì)最陡爬坡實(shí)驗(yàn)來快速逼近各顯著參數(shù)的最佳區(qū)域。根據(jù)小麥-小麥靜摩擦系數(shù)、小麥-小麥滾動摩擦系數(shù)和小麥-有機(jī)板靜摩擦系數(shù)3個因素的效應(yīng)大小設(shè)計(jì)相應(yīng)的步長，實(shí)驗(yàn)設(shè)計(jì)及結(jié)果如表5所示。在第二組時，休止角的相對誤差最小，因此選擇第二組為水平中心點(diǎn)，第一組為低水平，第三組為高水平進(jìn)行響應(yīng)面優(yōu)化實(shí)驗(yàn)。

表5 最陡爬坡實(shí)驗(yàn)設(shè)計(jì)及結(jié)果

2.3.3 響應(yīng)面優(yōu)化實(shí)驗(yàn)

根據(jù)Plackett-Burman和最陡爬坡實(shí)驗(yàn)結(jié)果，以小麥-小麥靜摩擦系數(shù)、小麥-小麥滾動摩擦系數(shù)和小麥-有機(jī)板靜摩擦系數(shù)為自變量，休止角θ為因變量建立三因素三水平實(shí)驗(yàn)，如表6所示。

表6 Box-Behnken實(shí)驗(yàn)設(shè)計(jì)及結(jié)果

通過Design-Expert求解得到休止角與3個顯著參數(shù)的二次多項(xiàng)式回歸模型為：

θ=26.51+4.13A+1.14B+3.17C-0.19B-1.18AC-1.05BC-1.32+0.37B2-4.57C2

式中：θ為休止角；A為小麥-小麥靜摩擦系數(shù)；B為小麥-小麥滾動摩擦系數(shù)；C為小麥-有機(jī)板靜摩擦系數(shù)。

回歸模型方差分析結(jié)果如表7。其中，回歸模型P<0.000 1(極其顯著)，失擬項(xiàng)P=0.52>0.05(不顯著)，說明該回歸模型可靠，擬合程度高。由A、B、C的P值可判斷出：三者對休止角的影響都極其顯著。此外，模型的決定系數(shù)R2=0.983 8，說明回歸方程相關(guān)性較好；校正決定系數(shù)adjR2=0.962 9，說明預(yù)測值與實(shí)際值具有較高的相關(guān)性；變異系數(shù)CV=3.72%，說明模型合理。此模型可用于對休止角的影響因素進(jìn)行分析。

表7 回歸模型方差分析結(jié)果

運(yùn)用Design-Expert中的Optimization功能，以實(shí)際休止角為目標(biāo)值，求解得到最優(yōu)參數(shù)為小麥-小麥靜摩擦系數(shù)為0.185、小麥-小麥滾動摩擦系數(shù)為0.046、小麥-有機(jī)板靜摩擦系數(shù)數(shù)為0.203。

3 驗(yàn)證實(shí)驗(yàn)

分別用構(gòu)建的4種小麥籽粒模型對標(biāo)定的離散元參數(shù)進(jìn)行仿真，并于實(shí)際實(shí)驗(yàn)進(jìn)行對比分析，得到的仿真結(jié)果結(jié)果及相對誤差如表8所示。從表8中可以看出，4種模型的仿真結(jié)果相近，與實(shí)際休止角的相對誤差小于1.30%，表明本文基于X射線斷層掃描方法構(gòu)建的模型、標(biāo)定的參數(shù)可用于離散元仿真。

表8 實(shí)驗(yàn)驗(yàn)證結(jié)果

4 結(jié)論

基于X射線斷層掃描技術(shù)，運(yùn)用圖像濾波、圖像分割等算法重構(gòu)了小麥籽粒的外輪廓模型；在EDEM中，分別填充4種不同數(shù)量的球(17、49、80、126球)來擬合小麥籽粒的真實(shí)形狀，構(gòu)建了更加精準(zhǔn)的小麥離散元模型。

基于X射線斷層掃描技術(shù)構(gòu)建的離散元仿真模型，通過休止角實(shí)驗(yàn)與離散元仿真實(shí)驗(yàn)結(jié)合，運(yùn)用Plackett-Burman實(shí)驗(yàn)、最陡爬坡實(shí)驗(yàn)、響應(yīng)面優(yōu)化實(shí)驗(yàn)，標(biāo)定了小麥離散元仿真所需要的參數(shù)；用構(gòu)建的4種小麥籽粒模型對標(biāo)定的離散元參數(shù)進(jìn)行仿真分析并與實(shí)際實(shí)驗(yàn)進(jìn)行驗(yàn)證，結(jié)果表明本文方法構(gòu)建的離散元模型及標(biāo)定的參數(shù)可用于離散元仿真。