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

李坤由，潘 樊，孫慧男

（鄭州中糧科研設(shè)計(jì)院有限公司，河南 鄭州 450001）

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

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

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

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

1 材料與方法

1.1 實(shí)驗(yàn)?zāi)Ｐ?/h3>

1.1.1 構(gòu)建模型

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

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

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

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

1.1.2 顆粒接觸模型

用干燥法測(cè)得該品種玉米的含水率為14.23%，由于其水含量較低，可忽略顆粒間的黏附力。假設(shè)顆粒堆積運(yùn)動(dòng)過程中，其應(yīng)力、位移、速度是通過顆粒與接觸體之間產(chǎn)生的微小交疊量值的不同來確定的?；诖思僭O(shè)，在EDEM軟件中建立了與相對(duì)速度相關(guān)的 Hertz-Mindlin（no slip）接觸力學(xué)模型來定義顆粒的接觸。該模型將顆粒在接觸點(diǎn)的碰撞力分解為法向力和切向力，如圖 3所示，切向和法向的接觸力簡(jiǎn)化成彈簧和阻尼器的并聯(lián)。

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

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

1.1.3 材料性能參數(shù)

玉米顆粒的泊松比按照Tao等[15]描述的結(jié)果作為參考值。依據(jù)標(biāo)準(zhǔn)ASAE S 368.4(2006)[16]，對(duì)單個(gè)顆粒進(jìn)行剪切實(shí)驗(yàn)分析，確定玉米的剪切模量，見表1。

表1 材料物理特性參數(shù)Table 1 Material physical properties parameters

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

表2 接觸參數(shù)Table 2 Contact parameter

1.2 實(shí)驗(yàn)方法

1.2.1 實(shí)驗(yàn)過程

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

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

圖4 實(shí)驗(yàn)流程圖Fig.4 Test flow chart

為減少實(shí)驗(yàn)的偶然性誤差，進(jìn)行多次重復(fù)實(shí)驗(yàn)。

1.2.2 模擬的合理性驗(yàn)證

將 1.1.3中列舉的材料性能參數(shù)作為基礎(chǔ)參數(shù)輸入到EDEM軟件中，進(jìn)行玉米顆粒的裝卸模擬，進(jìn)行五次平行重復(fù)實(shí)驗(yàn)，取其平均結(jié)果與實(shí)驗(yàn)平均結(jié)果做對(duì)比。

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

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

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

2 結(jié)果與分析

2.1 宏觀現(xiàn)象分析

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

增大顆粒間摩擦?xí)趦蓚€(gè)方面影響最終時(shí)間步長(zhǎng)。首先，它會(huì)產(chǎn)生較小的體積密度值，從而間接導(dǎo)致T的增加；此外，Balevicius[17]等發(fā)現(xiàn)顆粒間的滾動(dòng)摩擦?xí)绊?Beverroos[18]方程的C常數(shù)，從而影響卸料速度，即隨著摩擦的增加，卸料速率減小，卸料時(shí)間步長(zhǎng)增加。

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

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

2.2 顆粒力鏈的空間分布

2.2.1 切片參數(shù)

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

圖7 切片點(diǎn)位圖Fig.7 Location map of monitoring points

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

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

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

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

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

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

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

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

圖9為Z=4的力鏈分布情況。橫向?qū)Ρ炔煌瑫r(shí)間點(diǎn)的力鏈分布發(fā)現(xiàn)，從0.2 T到0.4 T不同摩擦系數(shù)下的玉米顆粒的強(qiáng)力連都呈增長(zhǎng)狀態(tài)，說明在Z=4高度處的玉米顆粒隨時(shí)間增加應(yīng)力在不斷聚集，未觀察到成拱塌陷現(xiàn)象。

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

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

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

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

3 結(jié)論

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

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

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

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