小麥籽粒振動篩分黏彈塑性接觸模型構建及其參數標定

劉凡一，張 艦，陳 軍

（西北農林科技大學機械與電子工程學院，楊凌 712100）

0 引 言

離散元法能夠實時地通過顆粒-顆粒/幾何體間的接觸來計算顆粒受力和運動，因此其在農業物料振動篩分中得到了廣泛應用[1-3]。顆粒接觸模型是離散元法的重要基礎，其準確性直接決定了顆粒所受力和力矩的大小，直接影響計算結果的準確性。

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

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

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

為更好地模擬麥粒振動篩分過程，本研究在現有離散元接觸模型基礎上，構建了一種黏彈塑性接觸模型，并結合真實試驗（加載-卸載、碰撞及旋轉鼓）對其參數進行標定，最后通過振動篩分試驗對標定結果進行驗證。

1 黏彈塑性接觸模型構建

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

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

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

其中

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

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

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

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

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

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

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

如果 Ft＞μFn，則

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

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

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

2 模型參數標定

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

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

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

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

2.1 單軸加載-卸載試驗

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

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

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

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

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

2.2 碰撞試驗

2.2.1 麥粒-麥粒

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

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

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

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

當麥粒屈服重疊量為7.63×10-6m時，改變麥粒-麥粒法向阻尼系數Cn-pp進行碰撞模擬，得到麥粒-麥粒恢復系數隨法向阻尼系數變化趨勢如圖5a所示。利用Origin軟件進行擬合，得到二者關系如下

以麥粒-麥粒恢復系數真實測量結果為目標值，求解二次方程可得麥粒-麥粒法向阻尼系數Cn-pp為190.68，此時模擬所得麥粒-麥粒恢復系數（0.44）與真實試驗值的相對誤差為4.35%。

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

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

2.2.2 麥粒-鋼板

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

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

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

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

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

2.3 旋轉鼓試驗

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

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

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

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

表1 旋轉鼓模擬試驗設計及其結果Table 1 Design and results of rotating drum simulation test

3 試驗驗證

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

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

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

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

4 結 論

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