基于PFC3D的玉米籽粒離散元參數(shù)研究

戴繼光

(沈陽農(nóng)業(yè)大學(xué)，遼寧 沈陽 110000)

我國農(nóng)業(yè)機(jī)械水平不斷提高，針對(duì)糧食物料特性也有了深入的研究。玉米是我國主要糧食作物之一，玉米堆可視為各向異性的離散體聚集而成的集合體，有限元、流體等分析方式不適用于處理離散物質(zhì)大量復(fù)雜行為信息，源于分子動(dòng)力學(xué)的離散元法是解決散體顆粒的內(nèi)部運(yùn)動(dòng)和受力情況的主要研究方法[1]。目前離散元已經(jīng)應(yīng)用于不同領(lǐng)域的多方面。農(nóng)業(yè)工程方面主要應(yīng)用于筒倉卸料、顆粒與粉體加工等[2]。許多學(xué)者將離散元法應(yīng)用于農(nóng)業(yè)生產(chǎn)、機(jī)械等方面用于分析物料在農(nóng)業(yè)生產(chǎn)過程的運(yùn)動(dòng)狀態(tài)[3-7]。物料的離散元參數(shù)是研究的一個(gè)重點(diǎn)[8-11]，深入研究玉米顆粒之間復(fù)雜信息需要得到更貼合實(shí)際的離散元參數(shù)，基于PFC3D軟件的玉米籽粒參數(shù)尚不完善，本文就解決這一問題做出研究。

于慶旭等采用EDEM軟件對(duì)三七種子進(jìn)行參數(shù)標(biāo)定[12]，模擬堆積試驗(yàn)得出三七種子離散元參數(shù)，堆積試驗(yàn)的靜止角反應(yīng)是糧堆摩擦，摩擦又分為內(nèi)摩擦和外摩擦。內(nèi)摩擦可以用內(nèi)摩擦角和靜止角衡量，內(nèi)摩擦角是用以計(jì)算倉儲(chǔ)設(shè)備重要力學(xué)參數(shù)[13]，靜止角試驗(yàn)進(jìn)行參數(shù)標(biāo)定也是離散元參數(shù)研究的一個(gè)重要方向[14-16]。直剪試驗(yàn)一直是研究土體抗剪強(qiáng)度、測定內(nèi)摩擦角和內(nèi)聚力等參數(shù)的有效途徑。唐福元等利用直剪儀對(duì)不同含水量的玉米的內(nèi)摩擦角進(jìn)行了試驗(yàn)測定得到準(zhǔn)確的內(nèi)摩擦角范圍[17]。室內(nèi)真實(shí)直剪試驗(yàn)只能獲得宏觀力學(xué)表現(xiàn)，史乃偉等采用PFC3D對(duì)沙土進(jìn)行了直剪試驗(yàn)?zāi)M以及力鏈分析[18]。離散元軟件PFC3D中不同的接觸模型下需要標(biāo)定的參數(shù)不同[19]。

本文基于直剪試驗(yàn)用PFC3D程序?qū)θS玉米顆粒進(jìn)行填充，著重研究了線性接觸模型之下在模擬直剪試驗(yàn)中對(duì)玉米對(duì)內(nèi)摩擦角影響較大的因子。以內(nèi)摩擦角為響應(yīng)值采用中心組合設(shè)計(jì)試驗(yàn)分析出其回歸方程，通過比對(duì)真實(shí)直剪試驗(yàn)的數(shù)據(jù)得出較為準(zhǔn)確的模擬參數(shù)。

1 材料與方法

1.1 直剪試驗(yàn)方法

玉米對(duì)的內(nèi)摩擦角采用直剪儀測定法，此方法是基于Coulomb理論的一種試驗(yàn)方法。直剪儀主要由剪切盒和伺服裝置組成，剪切盒分為上下兩部分，試驗(yàn)過程上盒固定并且施加一個(gè)垂直壓力，伺服裝置在下盒施加一個(gè)水平推力使得試驗(yàn)完全破壞。依據(jù)莫爾理論：

τ=c+σtgφ

(1)

式中，τ為剪應(yīng)力，kPa；σ為正應(yīng)力，kPa；φ為摩擦角，°；c為內(nèi)聚力，kPa。

假設(shè)糧食不具有內(nèi)聚力，得出內(nèi)摩擦角與應(yīng)力關(guān)系式：

φ=arctan(τ/σ)

(2)

1.2 室內(nèi)真實(shí)直剪試驗(yàn)

試驗(yàn)采用玉米品種為“遼單502”，含水量11.02%，密度1197kg·m-3。依據(jù)剪切盒與試料顆粒的比例關(guān)系選擇剪切盒凈空尺寸為100mm的正方體剪切盒[20]，剪切速率參考其它非粘性體快剪試驗(yàn)選為0.5mm·s-1,剪切過程壓力應(yīng)根據(jù)糧食30m深度的受力載荷為依據(jù)，取最大σ為240kPa，載荷等級(jí)應(yīng)該是50kPa、100kPa、150kPa、200kPa。利用土工合成材料綜合測定儀在4個(gè)不同壓力下重復(fù)測定5次得到的內(nèi)摩擦角的平均值見表1。

表1 室內(nèi)直剪試驗(yàn)結(jié)果

1.3 數(shù)值直剪試驗(yàn)

1.3.1 種粒模型建立

常見的玉米顆?？梢苑譃?種，體積較大的類馬齒形、體積中等的類錐形和類球形。隨機(jī)選取實(shí)體試驗(yàn)所用的玉米顆粒100粒，玉米粒下底>5mm視為類馬齒型，<5mm視為類錐型，體型較小且類似球體的視為類球型，按照上述標(biāo)準(zhǔn)分類并統(tǒng)計(jì)尺寸數(shù)據(jù)如表2。按照尺寸表進(jìn)行三維建模并導(dǎo)入PFC3D軟件進(jìn)行填充[21]。

表2 玉米籽粒參數(shù)

圖1 玉米籽粒填充模型

1.3.2 線性接觸模型

線性接觸模型的結(jié)構(gòu)如圖2所示，接觸力可分為線性部分F1和阻尼部分Fd。線性部分提供線彈性、摩擦行為，阻尼部分提供粘性行為。線性力通過具有恒定剛度的線性彈簧產(chǎn)生(ks、kn)。線性彈簧不能滿足張力，通過摩擦系數(shù)μ對(duì)剪切力施加庫倫準(zhǔn)則滿足滑移條件，阻尼力由阻尼器產(chǎn)生。

圖2 線性接觸模型示意圖

數(shù)值直剪試驗(yàn)剪切盒尺寸為100mm×100mm×100mm,在后續(xù)標(biāo)定試驗(yàn)中使豎直方向壓力為200kPa，經(jīng)循環(huán)后達(dá)到給定壓力狀態(tài)，試樣在穩(wěn)定壓力下剪切設(shè)定速率為0.06mm·s-1,剪切位移設(shè)定20mm,對(duì)整個(gè)過程的受力狀態(tài)進(jìn)行檢測。剪切過程如圖3所示。

圖3 數(shù)值直剪試驗(yàn)過程

1.3.3 試驗(yàn)設(shè)計(jì)方案

通過PB試驗(yàn)確定玉米籽粒參數(shù)中顯著因子，根據(jù)線性接觸后模型仿真需要的相關(guān)參數(shù)如表3所示。選取8個(gè)變量，每個(gè)參數(shù)選取高、低2個(gè)水平，以內(nèi)摩擦角為響應(yīng)，通過12組PB試驗(yàn)選擇3個(gè)顯著的參數(shù)。

表3 仿真參數(shù)設(shè)定

經(jīng)過PB試驗(yàn)選取的3個(gè)顯著參數(shù)變量，多次試驗(yàn)結(jié)果選取3個(gè)因素的中心點(diǎn)進(jìn)行三因素五水平的中心組合設(shè)計(jì)試驗(yàn)，對(duì)得到的試驗(yàn)結(jié)果進(jìn)行響應(yīng)面分析，得出回歸方程并繪出直觀的三維立體響應(yīng)面圖，通過計(jì)算預(yù)測出合適的3個(gè)參數(shù)。

為驗(yàn)證得到參數(shù)的有效性，將得到的參數(shù)再次進(jìn)行數(shù)值直剪試驗(yàn)，并且與室內(nèi)真實(shí)直剪試驗(yàn)結(jié)果比對(duì)兩者之間的誤差。

2 仿真試驗(yàn)結(jié)果

2.1 PB試驗(yàn)結(jié)果

根據(jù)PB試驗(yàn)的分析得出3個(gè)顯著的因子是有效模量阻尼系數(shù)和種間的摩擦系數(shù)，因此在后續(xù)試驗(yàn)中采用三因子中心組合設(shè)計(jì)試驗(yàn)。模擬試驗(yàn)結(jié)果如表4所示，由試驗(yàn)結(jié)果可知，第10組試驗(yàn)結(jié)果相較于室內(nèi)真實(shí)直剪試驗(yàn)的結(jié)果誤差最小，因此取第10組3顯著因子的值作為中心點(diǎn)，取玉米有效模量為4.09,阻尼系數(shù)為0.72,種間摩擦系數(shù)為0.55。后續(xù)中心組試驗(yàn)中其它非顯著因子也保持第10組PB試驗(yàn)中的數(shù)據(jù)。

圖4 PB試驗(yàn)帕累托圖

表4 Placket-Burman試驗(yàn)設(shè)計(jì)及結(jié)果

2.2 中心組合試驗(yàn)結(jié)果

依據(jù)PB試驗(yàn)得到的中心點(diǎn)進(jìn)行三因素五水平的中心組合設(shè)計(jì)試驗(yàn)。試驗(yàn)設(shè)計(jì)以及試驗(yàn)結(jié)果如表5所示。

根據(jù)試驗(yàn)結(jié)果，進(jìn)行回歸方程擬合，得出關(guān)于內(nèi)摩擦角和3因子之間的二次回歸方程：

(3)

利用Design-Expert軟件分析該模型的方差，結(jié)果如表5，得到該模型P=0.0003，說明該模型自變量與因變量相關(guān)性極其顯著，方程的擬合度極好，而擬失項(xiàng)P=0.0147，差異顯著，必須對(duì)該方程進(jìn)行手動(dòng)優(yōu)化。

表5 中心組合試驗(yàn)及結(jié)果

依據(jù)方差分析在原有二次回歸方程上增加高階項(xiàng)ABC、A2B、A2C后發(fā)現(xiàn)，該模型P<0.0001,該方程的擬失項(xiàng)P=0.3202，相關(guān)系數(shù)R2=0.9867校正決定系數(shù)R2Adj=0.9640，相較于優(yōu)化以前有較大提升，新擬合方程：

(4)

表6 優(yōu)化后方差分析

新擬合的模型A有效模量和B阻尼系數(shù)交互作用明顯。且B阻尼系數(shù)與C種間摩擦系數(shù)交互作用明顯。應(yīng)用Design Expert軟件繪制響應(yīng)曲面，圖5a顯示有效模量與阻尼系數(shù)的效應(yīng)面相對(duì)于阻尼系數(shù)有效模量的影響較大，而圖5b顯示阻尼系數(shù)與種間摩擦系數(shù)的效應(yīng)面類似于馬鞍面，兩者的交互作用明顯。

圖5 有效模量、阻尼系數(shù)與種間摩擦系數(shù)效應(yīng)面

內(nèi)摩擦角的目標(biāo)值為23.30°，應(yīng)用分析軟件對(duì)新模型預(yù)測得到最優(yōu)組合：玉米剪切模量為385kPa、阻尼系數(shù)為0.62、種間摩擦系數(shù)為0.051。

圖6 不同法向載荷下數(shù)值模擬試驗(yàn)的應(yīng)力變化圖

數(shù)值模擬試驗(yàn)和真實(shí)直剪試驗(yàn)的對(duì)比如表7所示，隨著法向載荷的變化，室內(nèi)真實(shí)直剪試驗(yàn)和數(shù)值模擬試驗(yàn)的結(jié)果誤差不大，均呈現(xiàn)隨著法向載荷增大內(nèi)摩擦角逐漸減小的趨勢，數(shù)值模擬直剪試驗(yàn)剪切應(yīng)力曲線趨勢相似。整體而言，試驗(yàn)結(jié)果充分證明響應(yīng)面法標(biāo)定參數(shù)的可行性，也證明了數(shù)值模擬試驗(yàn)的正確性。

表7 驗(yàn)證試驗(yàn)分析

3 總結(jié)

通過響應(yīng)面試驗(yàn)篩選出在線性接觸模型中對(duì)試驗(yàn)內(nèi)摩擦角影響顯著的3個(gè)因子分別為玉米彈性模量、粘性阻尼系數(shù)、玉米顆粒種間摩擦系數(shù)。

通過響應(yīng)面試驗(yàn)標(biāo)定得到的3個(gè)顯著因子參數(shù)為玉米剪切模量為385kPa、阻尼系數(shù)為0.62、種間摩擦系數(shù)為0.051，不顯著因子玉米與鋼板之間的摩擦系數(shù)為對(duì)比之間試驗(yàn)的結(jié)果采用響應(yīng)面法標(biāo)定玉米離散元參數(shù)可行，玉米顆粒之間應(yīng)不存在內(nèi)聚力，而仿真結(jié)果得出存在內(nèi)聚力的原因是玉米顆粒之間存在相互咬合自鎖行為。