基于滾動摩擦作用下玉米顆粒在模擬筒倉卸料中的力鏈分析

李坤由，潘 樊，孫慧男

（鄭州中糧科研設計院有限公司，河南 鄭州 450001）

糧食屬于散體顆粒，其摩擦特性要比工程材料之間或固體干摩擦復雜得多，特別是糧食顆粒的不規(guī)則性對其在筒倉卸料動態(tài)過程中所呈現的摩擦特性有很大影響，而顆粒的摩擦特性又直接或間接作用于顆粒的流動狀態(tài)，從而影響著糧倉的儲、卸安全。

在卸糧過程中，顆粒的流態(tài)及應力等散體動力學特征變化對筒倉安全具有重要意義。Weinhart[1]運用離散元法捕捉了 DEM 模擬中的顆粒應力；R. Koby?ka[2]通過模擬儲料卸糧過程中顆粒的速度、應力分布探究側壁壓力波動和摩擦特性的演化機制。Yung[3]通過實驗和模擬研究發(fā)現柔性球在剛性平面上運動時的滾動摩擦與滾動速度呈非線性關系。劉萬峰[4]，崔濤[5]，Ketterhagen[6]等設計了不同測定顆粒滾動摩擦的小型實驗并通過離散元法進行了驗證，為數值模擬提供了基礎參數。在顆粒流動現象的研究中，除了筒倉外形和導流體等外部因素，顆粒摩擦特性等內部因素的重要性也逐漸被重視。

Weizman[7]在進行數值模擬時提出將旋轉速度轉化為平移速度，確定了剛性球在柔性表面上滾動摩擦的變化系數。Silber[8]等研究確定顆粒堆中顆粒空間分布規(guī)律依賴顆粒摩擦系數。Goniva[9]發(fā)現在單噴口流化床仿真模擬中引入滾動摩擦后能夠極大地改善模擬結果。YongZeng[10]將卸糧速度場的波動與倉內力鏈周而復始的消失現象結合在一起。研究[11-14]將顆粒的宏觀現象與細觀變化對應起來，然而，即使在單一球模型情況下，筒倉中由于顆粒滾動摩擦耦合效應導致的顆粒力鏈分布變化和拱效應也沒有得到很好的解釋或研究。

基于以上分析，對不規(guī)則玉米顆粒在筒倉卸料動態(tài)過程的摩擦特性的研究具有重要意義。但流態(tài)作為可視化的宏觀變化現象，在實驗中很難對其進行定量分析，因此本文借助離散元建立模型研究卸糧動態(tài)過程中的細觀變量，搭建宏-細觀的橋梁，探究基于滾動摩擦作用下的不規(guī)則玉米顆粒在筒倉卸料動態(tài)模擬中力鏈演化機理，以期為筒倉的結構設計和優(yōu)化提供理論基礎。

1 材料與方法

1.1 實驗模型

1.1.1 構建模型

1.1.1.1 顆粒模型 本文模擬顆粒原型為我國種植面積較廣的德美亞一號玉米品種。樣品選取的玉米籽粒的三維幾何尺寸如圖1所示，上底寬I、下底寬L、高h。在模擬仿真中采用多球法近似玉米顆粒的形狀，共計六個球模型相互組合。

圖1 玉米顆粒外形尺寸圖Fig.1 Dimensions of corn grain

1.1.1.2 筒倉模型 實驗室的小型筒倉采用亞克力材料制作，見圖2a，得益于其良好的透明性，可以直接觀察到筒倉內顆粒在填充和排放過程中的運動。筒倉邊長為L=120 mm，倉高H=250 mm；其壁面與水平的夾角為60°，出料口邊長D=50 mm。按照實驗筒倉尺寸，采用3Dmax繪制同等尺寸的筒倉模型，如圖2b所示。

圖2 筒倉尺寸圖Fig.2 Diagram of silo size

1.1.2 顆粒接觸模型

用干燥法測得該品種玉米的含水率為14.23%，由于其水含量較低，可忽略顆粒間的黏附力。假設顆粒堆積運動過程中，其應力、位移、速度是通過顆粒與接觸體之間產生的微小交疊量值的不同來確定的。基于此假設，在EDEM軟件中建立了與相對速度相關的 Hertz-Mindlin（no slip）接觸力學模型來定義顆粒的接觸。該模型將顆粒在接觸點的碰撞力分解為法向力和切向力，如圖 3所示，切向和法向的接觸力簡化成彈簧和阻尼器的并聯。

圖3 顆粒接觸模型Fig.3 Particle contact model

其中，wi是角速度，rad/s；Tt是切向扭矩，N· m；Tr是滾動摩擦力矩，N· m；Rr是有效滾動接觸半徑，m；wrel是兩顆粒間的相對角速度，rad/s；rμ是滾動摩擦系數。

1.1.3 材料性能參數

玉米顆粒的泊松比按照Tao等[15]描述的結果作為參考值。依據標準ASAE S 368.4(2006)[16]，對單個顆粒進行剪切實驗分析，確定玉米的剪切模量，見表1。

表1 材料物理特性參數Table 1 Material physical properties parameters

通過借鑒Gonzalez[11]設計的測定玉米顆粒間以及顆粒與接觸材料間靜摩擦系數和恢復系數的裝置，測得玉米顆粒間及其與亞克力材料間的靜摩擦系數，見表2。

表2 接觸參數Table 2 Contact parameter

1.2 實驗方法

1.2.1 實驗過程

以黃色玉米籽粒為主材料，加入少量黑色玉米籽粒做分隔層，以便能夠清晰的觀察到卸糧過程中顆粒的流動狀態(tài)。采用中心裝料的方式裝載顆粒，在裝料過程中，由于顆粒的下落具有一定的速度，彼此之間會發(fā)生碰撞從而使得不同顏色顆粒發(fā)生輕微混合，故顏色界限并不完全水平。

裝料完成后，先將滿載的模型倉靜置2 h，然后打開數碼相機，抽走出料口處的擋板使顆粒自然下卸，記錄卸糧全過程，實驗流程如圖4所示。

圖4 實驗流程圖Fig.4 Test flow chart

為減少實驗的偶然性誤差，進行多次重復實驗。

1.2.2 模擬的合理性驗證

將 1.1.3中列舉的材料性能參數作為基礎參數輸入到EDEM軟件中，進行玉米顆粒的裝卸模擬，進行五次平行重復實驗，取其平均結果與實驗平均結果做對比。

從裝料方面看，筒倉內實驗裝料高度為245.2 mm，模擬裝料高度為250 mm；從卸糧時間來看，從抽出底板到顆粒完全流出漏斗，實驗的平均時長為6.3 s，模擬的平均時長為6.6 s，模擬誤差為4.7%。

對比實驗和模擬中顆粒的流動行為（圖5），可以觀察到，隨著卸糧進行，料床高度下降。在0.25T時刻，上部黑色顆粒層較為平直均勻，而下部呈現出明顯的V形，顆粒聚攏運動明顯。隨著時間的增加，筒倉上部黑色顆粒層彎曲明顯，在0.75T時，模擬及實驗的筒倉與料斗連接處都出現顆粒滯留現象，即出現死料區(qū)，整個卸糧過程中實驗與模擬的顆粒流動行為基本吻合，驗證了模擬的準確性。

圖5 實驗與模擬流態(tài)對比圖Fig.5 Comparison of experimental and simulated flow patterns

2 結果與分析

2.1 宏觀現象分析

倉壁與顆粒間的摩擦系數保持不變（為非零定值0.031 1），改變顆粒間的滾動摩擦系數μ，發(fā)現μ值不但會影響顆粒在筒倉中的流動形態(tài)，還對卸糧的時間步長有一定的影響。

增大顆粒間摩擦會在兩個方面影響最終時間步長。首先，它會產生較小的體積密度值，從而間接導致T的增加；此外，Balevicius[17]等發(fā)現顆粒間的滾動摩擦會影響 Beverroos[18]方程的C常數，從而影響卸料速度，即隨著摩擦的增加，卸料速率減小，卸料時間步長增加。

對比三種情況下的顆粒流動狀態(tài)（圖6），可以看出：摩擦系數越小，顆粒的運動越自由，不同顏色顆粒間的混合越劇烈，即在摩擦系數為 0狀態(tài)下，任意時刻，黑色顆粒“越界現象”更為明顯。且摩擦系數越小，黑色顆粒帶彎曲得越早，即顆粒的流動更早顯現出管狀流的特征，在2 s、3.5 s時刻，摩擦系數為0的第一層黑色顆粒帶最早開始出現彎曲，而摩擦系數越大，其第一層黑色顆粒帶彎曲越不明顯，進入 5 s后，顆粒基本都已進入管狀流狀態(tài)。故對于有漏斗的筒倉來說，顆粒間摩擦的減少會改變整體流和管狀流之間的極限，從而增加管狀流的面積。

圖6 三種滾動摩擦系數下顆粒流態(tài)對比Fig.6 Comparison of particle flow regimes under three rolling friction coefficients

2.2 顆粒力鏈的空間分布

2.2.1 切片參數

在筒倉倉體中心處建立 X-Y-Z方向上7-1-10的檢測網格，用以檢測筒倉內部各空間方位上的顆粒流動狀態(tài)及其細觀參數變化。網格具體劃分情況及編號見下圖7。

圖7 切片點位圖Fig.7 Location map of monitoring points

2.2.2 力鏈的空間分布及概率分析

鑒于筒倉上部顆粒會快速流失且下部靠近漏斗處顆粒速度較大，取Z=3、4、5區(qū)域的顆粒接觸力進行分析比較。

圖8左側為0時刻顆粒法向力/顆粒平均法向力，比值f在Z=3處水平面的空間分布情況，右側為f的概率分布情況。

圖8 力鏈分布圖（0時刻）Fig.8 Distribution map of force chain(T=0)

顆粒力鏈在空間分布上呈現出中間力鏈較強，四周力鏈弱的特點，形成強力鏈的顆粒構成顆粒體的主體骨架，可承受較大的荷載，形成弱力鏈的顆粒作為必要的填充，保證了顆粒體系的整體穩(wěn)定性。

在右側f的概率分布上，f ＞ 1的部分為強力鏈，f ＜ 1的部分為弱力鏈，f =1附近概率密度較大。這是由于顆粒體內部存在拱效應，在上部顆粒自重下顆粒體內部出現各向異性的力鏈網絡，該力鏈網絡作為荷載的傳遞路徑有強弱之分。對應圖中的力鏈變化趨勢，法向力越接近平均法向力，即f越接近1，其分布越廣泛，分布概率越大，在法向力與平均法向力相等時其分布概率達到峰值；而在f =1兩側的弱力鏈和強力鏈在顆粒體中的分布均較少，使得遠離f =1時力鏈的分布概率較小。

此外所取的上下三個區(qū)域的強力鏈存在互補的現象，其強力鏈與弱力鏈的占比也隨顆粒所在區(qū)域高度的不同發(fā)生變化，對應圖中 Z=3、4、5時強力鏈分別占 64%、55%、37%，即顆粒越靠近筒倉底部（Z=3），其強力鏈占比越高，顆粒擠壓力越大，形成拱的概率越高。

改變摩擦系數后，不同高度處的玉米顆粒在不同時刻的力鏈空間分布及概率分布情況如下圖9、10所示。

圖9為Z=4的力鏈分布情況。橫向對比不同時間點的力鏈分布發(fā)現，從0.2 T到0.4 T不同摩擦系數下的玉米顆粒的強力連都呈增長狀態(tài)，說明在Z=4高度處的玉米顆粒隨時間增加應力在不斷聚集，未觀察到成拱塌陷現象。

圖9 力鏈分布圖（Z=4）Fig.9 Distribution map of force chain (Z=4)

改變摩擦系數，玉米顆粒在筒倉中的力鏈空間分布在不同時間點有明顯差別。0.2 T時刻，摩擦系數為0.06時顆粒的強力鏈分布最少，通過概率統(tǒng)計其中 f≥1的強力鏈占比25%，弱力鏈占據大部分剖面；隨時間增加到0.4 T，應力聚集，強力鏈占比達到50%，強力鏈增長速度達到100%。摩擦系數為0時，0.2 T時刻強力鏈占比為38%，較μ=0.06時增長了52%，說明在這一高度降低摩擦系數玉米顆粒進行應力聚集的時間提前；隨著時間的增加，摩擦系數為0的玉米顆粒應力在不斷聚集，在0.4 T時達到53%，但其增長速度卻小于μ=0.06的情況，應力增長比較平穩(wěn)，不易出現成拱現象。摩擦系數為0.102時，0.2 T時刻強力鏈占比為48%，較μ=0.06時增長了92%；隨著時間的增加，強力鏈占比仍在不斷增大，在0.4 T時達到78%，遠高于其他兩種摩擦情況，但其增長速率上仍小于μ=0.06，說明在這一高度增大滾動摩擦系數玉米顆粒進行應力聚集的現象明顯，但應力增長速率相對并不高。

圖10為Z=3高度的力鏈分布情況，該高度位于倉身靠下 1/3處，由于更接近漏斗口，在前期卸料時流動截面收縮，顆粒迅速集聚，導致顆粒間應力增加，即在0.2 T時刻Z=3高度處的強力鏈占比都比同條件下Z=4高度的占比高。其中以μ=0.06條件下的力鏈變化最為顯著，強力鏈由25%增長至 92%，顆粒迅速集聚，應力激增，隨著時間的增加，強力鏈占比減小，出現應力消散，形成了起拱-塌陷效應。削弱顆粒間滾動摩擦μ=0，顆粒的強力鏈占比由Z=4高度的38%到Z=3高度的 49%，應力有所增長，但其增長幅度遠小于μ=0.06條件下的應力變化；隨著卸料的進行，強力鏈仍在小幅度增長，未出現應力消散，說明減小滾動摩擦玉米顆粒在Z=3到Z=4高度保持較穩(wěn)定的暢通卸料，未出現起拱的應力激增，以及拱塌陷的應力衰減。這一應力變化現象與該摩擦條件下卸料時長最短相吻合。增大顆粒間滾動摩擦μ=0.102，在0.2 T時刻，強力鏈的占比隨其下落高度的變化基本保持不變，隨著卸料的進行，0.4 T時刻強力鏈占比為31%，相比0.2 T時刻的49%減少了37%，出現應力消散現象，但與其他兩種摩擦條件不同的是，該時刻Z=3高度的強力鏈相比Z=4高度的反而有所減小，由31%減少至28%，說明增大顆粒間滾動摩擦不但會增加供效應，且堆積成拱高度距離漏斗口更高。

圖10 力鏈分布圖（Z=3）Fig.10 Distribution map of force chain (Z=3)

3 結論

本文利用EDEM建立了筒倉模型和玉米仿真顆粒模型，模擬裝卸糧實驗，并以物理卸糧實驗驗證了模型及模擬的準確性。在此基礎上對模擬數據進行分析，探究了不同滾動摩擦系數條件下，玉米顆粒卸料時的力鏈分布規(guī)律，分析了筒倉內起拱-塌陷現象的應力變化。模擬得出以下結論：

（1）改變顆粒間的滾動摩擦系數，會影響顆粒在筒倉中的流動形態(tài)。顆粒間滾動摩擦系數越小，顆粒的流動由整體流轉變?yōu)楣軤盍鞯臅r間越早。對于有漏斗的筒倉來說，顆粒間摩擦的減少會改變整體流和管狀流之間的極限，從而增加產生管狀流的面積。

（2）0 T時刻筒倉中下部三個區(qū)域的強力鏈存在互補的現象，其強力鏈與弱力鏈的占比也隨顆粒所在區(qū)域高度的不同發(fā)生變化，即顆粒越靠近筒倉底部，其強力鏈占比越高，顆粒擠壓力越大，形成拱的概率越高。

（3）卸料初期，摩擦系數為0.06時顆粒的強力鏈分布最少，但隨時間增加，玉米顆粒迅速集聚，應力增幅最大，后期，隨著顆粒高度的減小，出現應力消散，呈現起拱-塌陷效應；降低摩擦系數，玉米顆粒進行應力聚集的時間提前，應力增長比較平穩(wěn)，未出現起拱的應力突增，以及拱塌陷的應力衰減；增大滾動摩擦系數，玉米顆粒進行應力集聚的現象明顯，滾動摩擦不但會增加供效應，且出現成拱高度距離漏斗口更高。