小麥籽粒振動篩分黏彈塑性接觸模型構(gòu)建及其參數(shù)標(biāo)定

劉凡一，張 艦，陳 軍

（西北農(nóng)林科技大學(xué)機(jī)械與電子工程學(xué)院，楊凌 712100）

0 引 言

離散元法能夠?qū)崟r地通過顆粒-顆粒/幾何體間的接觸來計算顆粒受力和運(yùn)動，因此其在農(nóng)業(yè)物料振動篩分中得到了廣泛應(yīng)用[1-3]。顆粒接觸模型是離散元法的重要基礎(chǔ)，其準(zhǔn)確性直接決定了顆粒所受力和力矩的大小，直接影響計算結(jié)果的準(zhǔn)確性。

離散元接觸模型主要包括兩類：黏彈性接觸模型和滯回接觸模型。其中根據(jù)顆粒加載過程中的力-位移關(guān)系，黏彈性接觸模型又可分為線性黏彈性和非線性黏彈性接觸模型。1979年，Cundall和Strack最早提出linear spring模型[4]，該模型具有計算量小，計算時步大的優(yōu)點。然而，linear spring模型采用線性彈簧和阻尼，其力-位移關(guān)系是完全線性的，這與大多數(shù)生物質(zhì)材料非線性特性不符；此外，該模型在計算法向彈簧剛度時，需要指定顆粒特征速度，這增加了模型的使用難度。為解決力-位移關(guān)系完全線性的問題，學(xué)者們根據(jù)Hertz非線性接觸理論先后提出了多種非線性黏彈性接觸模型，其中常用的包括：Hertz-Mindlin(no slip)模型[5]、Kuwabara and Kono 模型[6]。

為模擬顆粒材料的塑性變形，Walton等[7-11]先后提出了多種滯回接觸模型。根據(jù)力-位移關(guān)系的不同，滯回接觸模型可分為線性滯回模型和非線性滯回模型。滯回接觸模型都有一個共同的特點，即加載和卸載階段彈簧剛度不同，以此表征接觸過程中的塑性變形及能量損失。

關(guān)于現(xiàn)有接觸模型在農(nóng)業(yè)物料顆粒動力學(xué)模擬中的計算精度，Wojtkowski等[12]通過不同含水率的油菜-鋼板碰撞研究了黏彈性接觸模型和滯回接觸模型的計算精度，結(jié)果表明：對于低含水率的顆粒，滯回接觸模型計算精度高于黏彈性接觸模型；但對于高含水率的顆粒，結(jié)果相反。Johnstone在應(yīng)用Hertz-Mindlin(no slip)模型模擬谷物顆粒固結(jié)試驗時發(fā)現(xiàn)：由于該接觸模型缺少塑性，導(dǎo)致卸載過程中發(fā)生過量的彈性回彈，同時這也使得再加載-卸載過程中不再有形變發(fā)生[13]。與大部分工業(yè)生產(chǎn)中的顆粒材料不同，農(nóng)業(yè)固體物料流變學(xué)特性研究表明農(nóng)業(yè)物料顆粒往往同時具有 3種基本性質(zhì)：彈性、黏性以及塑性[14]?，F(xiàn)階段，有關(guān)黏彈塑性材料接觸理論已進(jìn)行了大量研究[15-17]，但對農(nóng)業(yè)物料黏彈塑性接觸模型的研究還較少。盧立新等[18]基于果實等速壓縮變形特征,提出了由3個非線性彈簧元件、1個黏性元件以及1個滑塊元件組成的黏彈塑性模型，并通過試驗對所建立的模型進(jìn)行了驗證。鄭曉等[19-20]在物料流變學(xué)試驗基礎(chǔ)上，采用理論模型與經(jīng)驗?zāi)Ｐ徒Y(jié)合的方法，先后建立了油菜籽、芝麻及花生的非線性黏彈塑性本構(gòu)模型。

為更好地模擬麥粒振動篩分過程，本研究在現(xiàn)有離散元接觸模型基礎(chǔ)上，構(gòu)建了一種黏彈塑性接觸模型，并結(jié)合真實試驗（加載-卸載、碰撞及旋轉(zhuǎn)鼓）對其參數(shù)進(jìn)行標(biāo)定，最后通過振動篩分試驗對標(biāo)定結(jié)果進(jìn)行驗證。

1 黏彈塑性接觸模型構(gòu)建

該黏彈塑性接觸模型通過 EDEM 2.7軟件應(yīng)用程序編程接口（application programming interface, API）實現(xiàn)，該接口允許用戶通過C++語言構(gòu)建用戶自定義接觸模型、自定義顆粒工廠以及自定義耦合接口。

法向方向，該接觸模型通過將Kuwabara and Kono非線性黏彈性接觸模型中的黏性耗散項引入Thornton滯回接觸模型[9]進(jìn)行構(gòu)建。與 Thornton接觸模型相似，該接觸模型法向分量包括4個階段。

第一階段為黏彈性加載/卸載階段，該階段顆粒接觸重疊量小于屈服重疊量δy，其法向接觸力Fn計算式為

其中

式中E*為等效彈性模量，Pa；R*為等效接觸半徑，m；δn為顆粒-顆粒/幾何體法向接觸重疊量，m；Cn為顆粒-顆粒/幾何體法向阻尼系數(shù)；為顆粒-顆粒/幾何體法向相對速度，m/s；Ei、Ej為顆粒/幾何體的彈性模量（i、j均為顆粒/幾何體編號），Pa；νi、νj為顆粒/幾何體的泊松比；Ri、Rj為顆粒/幾何體的接觸半徑，m。

第二階段為黏塑性加載階段，其法向接觸力Fn計算式為

第三階段為黏彈性卸載/重加載階段，其法向接觸力Fn計算式為

式中Fnmax為最大法向接觸力，N；Fny為屈服力，其值等于顆粒發(fā)生屈服時的法向接觸力，N；δnmax為最大法向接觸重疊量，m。

第四階段為顆粒塑性變形恢復(fù)階段，顆粒間法向接觸力Fn始終為0，直至顆粒發(fā)生完全分離。

切向方向采用簡化Thornton切向接觸模型[21]。該接觸模型基于Mindlin and Deresiewicz非滑移接觸理論，計算切向力時在每一時步中對上一時步切向力進(jìn)行修正，第n時步切向接觸力計算式如式（6）～（8）所示。如果 Ft≤μFn，則

式中ΔFn為法向接觸力增量，N；kt為切向彈簧剛度，N/m；Δδt為切向重疊量增量，m；a為接觸面半徑，m；G*為等效剪切模量，Pa；Gi、Gj為顆粒/幾何體的剪切模量，Pa。

如果 Ft＞μFn，則

式中μ為顆粒-顆粒/幾何體靜摩擦系數(shù)。

滾動摩擦力矩計算同 Hertz-Mindlin(no slip)接觸模型，通過接觸面上的力矩Ti計算得到，即

式中μr為顆粒-顆粒/幾何體滾動摩擦系數(shù)；ωi為接觸點處顆粒的單位角速度矢量，rad/s。

2 模型參數(shù)標(biāo)定

標(biāo)定試驗中所用麥粒由西北農(nóng)林科技大學(xué)小麥育種中心提供，品種號為西農(nóng) 223，含水率為 10.05%，密度為1 350.30 kg/m3。麥粒為三軸不等顆粒，隨機(jī)選取150顆利用數(shù)顯游標(biāo)卡尺（精度：0.01 mm）對其三軸尺寸進(jìn)行測量，得到平均長、寬、厚分別為6.21、3.16和2.93 mm。麥粒彈性模量通過 ASAE S368.4[22]標(biāo)準(zhǔn)測量得到，其值為1.93×109Pa，泊松比取0.40[23-25]。研究中幾何體材料為鋼，其密度、剪切模量和泊松比分別為7 800 kg/m3、7×1010Pa和 0.30。

大量研究表明，精確的顆粒模型并不能顯著提高模擬精度，通常只需通過少量的單元球建立粗略的非球形模型即可得到較好的模擬精度[26-28]。因此，根據(jù)麥粒三軸尺寸，采用如圖4所示的19球填充顆粒模型（長×寬×厚：6.21 mm×3.16 mm×2.93 mm）。在旋轉(zhuǎn)鼓及振動篩分模擬中，基于真實麥粒長軸尺寸分布，確定仿真中麥粒按照標(biāo)準(zhǔn)正態(tài)分布生成（平均值為1，標(biāo)準(zhǔn)差為0.06）；為避免生成過小及過大的顆粒，根據(jù)實際測量得到的最小、最大顆粒尺寸，將顆粒半徑限制在 0.65～1.21倍的初始半徑之間。在所有模擬中，仿真時步取1×10-7s，網(wǎng)格尺寸取最小球形單元半徑的3倍。

據(jù)上述黏彈塑性接觸模型中法向接觸力計算公式，可知待確定參數(shù)包括麥粒屈服重疊量 δy以及麥粒-麥粒/鋼板法向阻尼系數(shù)Cn。單軸加載-卸載試驗中，加載板速率較小，整個過程可看作是準(zhǔn)靜態(tài)的，可以忽略法向阻尼系數(shù)Cn的影響（模擬中阻尼系數(shù)為0），因此采用單軸加載-卸載試驗即可對麥粒屈服重疊量δy進(jìn)行標(biāo)定?；謴?fù)系數(shù)反應(yīng)了物體碰撞過程中的能量損耗，故選用碰撞試驗中所測恢復(fù)系數(shù)去標(biāo)定麥粒-麥粒/鋼板法向阻尼系數(shù)。

對于切向接觸力和滾動摩擦力矩，待確定參數(shù)包括麥粒-麥粒/鋼板靜摩擦系數(shù)及滾動摩擦系數(shù)。顆粒材料在旋轉(zhuǎn)鼓試驗中的動態(tài)休止角能夠很好地反映顆粒的流動特性，常被用于顆粒離散元參數(shù)的標(biāo)定[13]，此外麥粒在振動篩分過程中是動態(tài)的，因此選用旋轉(zhuǎn)鼓試驗對其摩擦系數(shù)進(jìn)行標(biāo)定。預(yù)試驗結(jié)果表明，對于本文中所建立的模型，滾動摩擦系數(shù)對動態(tài)休止角影響較小，因此本文中麥粒-麥粒/鋼板滾動摩擦系數(shù)均取 0.05，只對麥粒-麥粒/鋼板靜摩擦系數(shù)進(jìn)行標(biāo)定。

2.1 單軸加載-卸載試驗

隨機(jī)選取10顆麥粒，采用DDL10型電子萬能試驗機(jī)（長春機(jī)械科學(xué)研究院有限公司，長春）進(jìn)行單軸加載-卸載試驗（加載速度為1 mm/min），得到麥粒平均力-位移曲線如圖1所示。由力-位移曲線可以看出，即使在位移量較小的情況下（0.10 mm），麥粒仍然表現(xiàn)出顯著的塑性行為。

圖1 麥粒單軸加載-卸載試驗中力-位移曲線Fig. 1 Force-displacement curves in uniaxial loading-unloading test of wheat

利用EDEM軟件模擬麥粒單軸加載-卸載試驗，得到各采樣點處模擬-試驗接觸力誤差平方和D隨屈服重疊量δy的變化規(guī)律如圖2所示。通過Origin 8.0軟件擬合，發(fā)現(xiàn)二者關(guān)系可以通過方程（10）表示。

圖2 接觸力誤差平方和隨麥粒屈服重疊量變化趨勢Fig. 2 Changes of sum of squared errors of contact force with yield overlap of wheat

求解方程(10)可得，當(dāng)麥粒屈服重疊量為7.63×10-6m時，模擬-試驗接觸力誤差平方和最?。―=60.83 N），此時模擬所得麥粒力-位移曲線如圖1所示，從中可以看出模擬值與真實試驗值非常接近，說明麥粒屈服重疊量標(biāo)定結(jié)果良好。

2.2 碰撞試驗

2.2.1 麥粒-麥粒

麥粒-麥?；謴?fù)系數(shù)采用 Gonzálezmontellano等[29]所用雙擺裝置進(jìn)行測量，測量原理如圖3。測量時，為防止釋放時麥粒具有一定的初速度，顆粒 1由真空泵吸引，靜止在高度為h0處。關(guān)閉真空泵，顆粒1在重力作用下開始向左擺動，至最低點處（高度為0）與顆粒2發(fā)生碰撞。碰撞后，顆粒1、2均向左擺動，擺動高度分別為h1和h2，麥粒-麥?；謴?fù)系數(shù)epp可通過式（11）計算。

整個測量過程通過奧林巴斯I-speed3高速攝影機(jī)（奧林巴斯（中國）有限公司，廣州）以500 fps記錄。測量時，選取5組（10顆）麥粒，每組試驗重復(fù)20次，得到麥粒-麥粒恢復(fù)系數(shù)為0.46±0.09。

圖3 麥粒-麥粒雙擺碰撞示意圖Fig. 3 of double pendulum test of wheat-wheat collision

離散元模擬中，通過EDEM軟件自定義顆粒工廠接口（API），將麥粒-麥粒雙擺碰撞簡化為無重力水平碰撞。如圖4所示，碰撞前，顆粒1、2水平速度分別為v1、v2；碰撞后，兩顆粒的速度分別為v3、v4。根據(jù)恢復(fù)系數(shù)的定義，麥粒-麥?；謴?fù)系數(shù)epp計算如式（12）所示。

當(dāng)麥粒屈服重疊量為7.63×10-6m時，改變麥粒-麥粒法向阻尼系數(shù)Cn-pp進(jìn)行碰撞模擬，得到麥粒-麥?；謴?fù)系數(shù)隨法向阻尼系數(shù)變化趨勢如圖5a所示。利用Origin軟件進(jìn)行擬合，得到二者關(guān)系如下

以麥粒-麥?；謴?fù)系數(shù)真實測量結(jié)果為目標(biāo)值，求解二次方程可得麥粒-麥粒法向阻尼系數(shù)Cn-pp為190.68，此時模擬所得麥粒-麥?；謴?fù)系數(shù)（0.44）與真實試驗值的相對誤差為4.35%。

圖4 麥粒-麥粒雙擺試驗離散元模擬Fig. 4 DEM simulation of wheat-wheat double pendulum test

圖5 恢復(fù)系數(shù)隨法向阻尼系數(shù)變化趨勢Fig. 5 Changes of coefficient of restitution with normal damping coefficient

2.2.2 麥粒-鋼板

麥粒-鋼板恢復(fù)系數(shù)采用類似黃小毛等[30]所用傾斜平面法進(jìn)行測量，測量原理如圖6a所示。測量時，同樣采用真空泵通過吸頭將麥粒吸住，此時麥粒相對碰撞點 O的高度為H0。關(guān)閉真空泵，麥粒與傾斜板發(fā)生碰撞后落在落料板上，落料板與碰撞點O的相對高度為H1，麥粒水平位移為S1；改變落料板與碰撞點O的相對高度（H2），得到麥粒水平位移為 S2，此時麥粒與傾斜板間的恢復(fù)系數(shù)epw可通過式（14）計算。

式中von、vn分別為碰撞前、后法向速度分量，m/s；vx、vy分別為碰撞后水平、垂直速度分量，m/s；v0為碰撞前垂直速度分量，m/s。

圖6 麥粒-傾斜面碰撞試驗Fig. 6 Wheat-inclined plane collision test

為防止麥粒與落料板接觸后發(fā)生回彈，在落料板上涂上一層黃油以“捕獲”麥粒。同時，麥粒與傾斜板碰撞后極易產(chǎn)生旋轉(zhuǎn)，導(dǎo)致被測物同時具有X和Y方向位移，這與恢復(fù)系數(shù)計算時顆粒只具有 X方向的位移假設(shè)不符。因此，為減小麥粒 Y方向位移導(dǎo)致的測量誤差，黃油在Y方向上的寬度為50 mm，即麥粒Y向位移小于25 mm時，麥粒將被落料板“捕獲”，如圖6b所示。測量時，隨機(jī)選取50顆麥粒進(jìn)行測量，得到麥粒-鋼板恢復(fù)系數(shù)為0.48。通過EDEM軟件模擬顆粒-傾斜面碰撞過程，改變麥粒-鋼板法向阻尼系數(shù)Cn-pw，求得不同法向阻尼系數(shù)下的麥粒-鋼板恢復(fù)系數(shù)，如圖5b所示。擬合發(fā)現(xiàn)二者關(guān)系可以通過二次多項式進(jìn)行表示。

以麥粒-鋼板恢復(fù)系數(shù)真實測量結(jié)果為目標(biāo)值，求解所建立二次方程可得麥粒-鋼板法向阻尼系數(shù)為306.65。采用標(biāo)定參數(shù)模擬麥粒-傾斜面碰撞，得到恢復(fù)系數(shù)模擬值（0.479 02）與真實測量值（0.484 35）差距極小。

2.3 旋轉(zhuǎn)鼓試驗

旋轉(zhuǎn)鼓試驗設(shè)置與Johnstone研究中的設(shè)置相似[13]。其中，鋼質(zhì)旋轉(zhuǎn)鼓內(nèi)徑和厚度分別為150和30 mm，為觀察記錄旋轉(zhuǎn)鼓內(nèi)顆粒流動，旋轉(zhuǎn)鼓一側(cè)采用有機(jī)玻璃擋板，如圖7a所示。為獲得平整的顆粒自由流動表面，旋轉(zhuǎn)鼓的轉(zhuǎn)速設(shè)為7 r/min，顆粒填充率為40%。麥粒在旋轉(zhuǎn)鼓中的流動圖像通過數(shù)碼相機(jī)獲?。扉T1/1 600 s），隨機(jī)選取5張圖片進(jìn)行動態(tài)休止角測量。

圖7 旋轉(zhuǎn)鼓試驗及模擬對比Fig.7 Comparison of rotating drum in lab test and simulation

為減小人為因素導(dǎo)致的測量誤差，采用計算機(jī)圖像技術(shù)進(jìn)行動態(tài)休止角測量。測量時，首先采用Image J軟件截取圓形旋轉(zhuǎn)鼓圖像；然后在MATLAB中采用K均值聚類的方法提取圖像邊界點；最后通過線性擬合得到擬合邊界，該直線斜率的反正切值即為麥粒動態(tài)休止角，本研究中測得麥粒動態(tài)休止角為（33.87±0.40）°。

動態(tài)旋轉(zhuǎn)鼓模擬試驗設(shè)置與真實試驗一致，仿真總時間為5 s，從4.2 s開始，每隔0.2 s截取顆粒流動圖像。結(jié)合國內(nèi)外離散元模擬中麥粒-麥粒/鋼靜摩擦系數(shù)的取值[13,25,31-33]，確定了參數(shù)變化范圍，具體參數(shù)組合及模擬結(jié)果如表1所示。從表1中可知，當(dāng)麥粒-麥粒/鋼板靜摩擦系數(shù)分別為 0.40和 0.44時，旋轉(zhuǎn)鼓動態(tài)休止角為（34.79±0.37）°，最接近真實試驗值（如圖7b所示）。

表1 旋轉(zhuǎn)鼓模擬試驗設(shè)計及其結(jié)果Table 1 Design and results of rotating drum simulation test

3 試驗驗證

為驗證所建立的接觸模型及標(biāo)定的模型參數(shù)能否用于麥粒振動篩分試驗，選用魚鱗篩進(jìn)行篩分試驗，如圖8所示。振動篩分試驗中，魚鱗篩開角、篩面傾角及振動方向角分別為15°、2°和15°，振幅為30.71 mm，頻率為4 Hz。研究中篩分試驗重復(fù)3次，每次試驗中麥?？傎|(zhì)量為3.90 kg。模擬中考慮到篩面在寬度方向的完全對稱及計算量[3]，寬度方向采用周期性邊界條件選取100 mm進(jìn)行建模，仿真中共生成31 706顆麥粒。

圖8 振動篩分試驗及其離散元模擬Fig.8 Vibratory screening test and its DEM simulation

為比較振動篩分試驗及其離散元模擬結(jié)果，在篩面下方設(shè)置接料盒，并在長度方向?qū)⑵浞譃?7個區(qū)間。待篩分完成后，統(tǒng)計試驗及模擬所得各區(qū)間內(nèi)顆粒質(zhì)量分?jǐn)?shù)，結(jié)果如圖9所示。由圖9可知，試驗及模擬結(jié)果具有極為相似的變化趨勢，篩下物分布均主要存在于前 3個區(qū)間，最大誤差出現(xiàn)在第 2區(qū)間（8.97%）。為評價試驗與模擬結(jié)果有無顯著性差異，采用SPSS軟件對2組數(shù)據(jù)進(jìn)行配對T檢驗，得到其P值為0.99，遠(yuǎn)大于0.05，說明2組數(shù)據(jù)無顯著性差異的概率為99%。同時，2組數(shù)據(jù)的相關(guān)系數(shù)為0.98，說明兩組數(shù)據(jù)具有極高的相關(guān)性。綜合分析，兩組數(shù)據(jù)無顯著性差異且高度相關(guān)，這表明所建立的黏彈塑性接觸模型及標(biāo)定的參數(shù)能夠很好地模擬麥粒振動篩分過程。

圖9 振動篩分試驗及模擬中篩下物分布對比Fig.9 Comparison of distribution of wheat in vibratory screening test and its simulation

4 結(jié) 論

本文以麥粒為研究對象，在現(xiàn)有離散元接觸模型基礎(chǔ)上，建立了一種黏彈塑性接觸模型。利用單軸加載-卸載試驗、碰撞試驗及旋轉(zhuǎn)鼓試驗對麥粒接觸模型參數(shù)進(jìn)行標(biāo)定，得到麥粒屈服重疊量為7.63×10-6m、麥粒-麥粒/鋼板法向阻尼系數(shù)分別為 190.68和 306.65、麥粒-麥粒/鋼板靜摩擦系數(shù)的最佳組合為0.40和0.44。麥粒振動篩分驗證試驗表明所建立的接觸模型及標(biāo)定的模型參數(shù)能夠很好地模擬麥粒振動篩分過程，模擬與試驗所得篩下麥粒質(zhì)量分?jǐn)?shù)最大誤差為8.97%。研究結(jié)果可為其他農(nóng)業(yè)物料黏彈塑性模型構(gòu)建及其參數(shù)標(biāo)定提供參考。