基于EDEM的糧食立筒倉卸料動態側壓力研究

劉強, 張茜茜, 劉濤, 丁永剛, 許啟鏗, 孫啟帥

(1.河南工業大學土木工程學院, 鄭州 450001; 2.河南省糧油倉儲建筑與安全重點實驗室, 鄭州 450001;3.機械工業第六設計研究院有限公司, 鄭州 450001)

立筒倉以其占地面積小、易于機械化、自動化作業、流通費用低、造價低等特點,被廣泛應用于糧食儲備。但貯料[1]與倉體存在著復雜的相互作用,貯料會受到倉壁的約束作用,在卸料過程中又會對倉壁產生一定的反作用力,發生超壓現象從而加大倉體破壞的風險,影響糧食的供給保障。

隨著筒倉需求量越來越大,相關的研究也在不斷完善。Wang等[2]運用細砂和玻璃微珠兩種散體材料,對比研究筒倉卸料時試驗和模擬的流動規律及法向壁壓分布情況,得出試驗和模擬的結果基本一致。Wójcik等[3]通過多次裝卸料試驗,得出在整體流狀態時,最大水平壓力和豎向壓力位于倉壁底部范圍內。Fullard[4]采用PIV對筒倉的重力驅動卸料試驗進行錄像分析,表明了在一個比例時間內,遠離排放孔,流動行為是相同的。原方等[5]對比分析普通筒倉安裝減壓管時筒倉單側壁卸料的動態側壓力,得出減壓管能夠減小卸料時儲料對倉壁的沖擊,從而提高卸料時筒倉的穩定性。王振清等[6]進行了50%、80%和100%倉容3 種狀態下的側壓力試驗,發現管狀流動的出現位置與初始的儲糧倉容相關,結果又與規范值進行對比,可為糧食筒倉設計提供參考。徐志軍等[7]分析顆粒在普通筒倉雙側壁卸料和帶流槽側壁卸料過程中的力學行為,驗證了溜槽的減壓機理。程遠浩[8]將曲線漏斗筒倉與錐形漏斗筒倉對比,曲線漏斗改善流態,從而降低了筒倉動態側壓力。石鑫等[9]建立的改進顆粒模型與傳統單元模型比較,發現其倉壁側壓力峰值顯著減小,力鏈分叉現象明顯。余漢華等[10]采用3種機器學習方法預測筒倉動態側壓力,同時進行隨機抽樣得到筒倉動態側壓力的概率分布,為筒倉結構的可靠度提供了理論基礎。李坤由等[11]通過EDEM(engineering discrete element method)離散元軟件,對比分析不同摩擦情況下力鏈的細觀參數,得出顆粒間摩擦系數影響卸料時間,增大顆粒間滾動摩擦會增加拱效應。

綜上可知,前人研究中針對卸料過程中形成的不同流動狀態沒有具體分析其成因;在靜態、動態側壓力等方面的研究成果雖有很多,但是對于超壓系數沒有固定的取值,且最大值范圍也并不確定;同時卸料過程中形成的動態拱具有連續形成和倒塌的特性,會對研究帶來較大的困難,有關動態拱的形成機理仍缺少相應解釋。為此,應用離散元軟件EDEM(engineering discrete element method),對比分析糧食立筒倉卸料時靜態、動態側壓力和超壓系數差異規律,探究筒倉卸料時成拱機理與倉壁側壓力的關系,為立筒倉動態側壓力提供一種新的研究途徑。

1 離散元模型建立

1.1 立筒倉模型

以某糧庫儲藏大豆的柱承式鋼筋混凝土立筒倉為研究對象,在本團隊已開展的試驗研究[12]基礎上,通過EDEM(engineering discrete element method)離散元軟件建立糧食立筒倉模型,采用有機玻璃材料屬性,模型總高度為1.62 m,其中倉壁高度1.28 m,內徑0.48 m,外徑0.50 m。結構模型物理參數泊松比為0.33,密度為1 180 kg/m3,彈性模量為3.05×103MPa,結構模型整體與生成貯料區域重力加速度為9.8 m/s2,糧食立筒倉模型如圖1所示。

圖1 立筒倉模型Fig.1 Silo model

1.2 貯料顆粒模型

所采用的貯料為陶粒,其為球形顆粒形態,在EDEM(Engineering Discrete Element Method)軟件前處理模塊的Bulk Material中創建貯料顆粒,選取Spheres,設置半徑為0.001 5 m,將顆粒的粒徑分布設置為正態分布,平均值設為1即表示所有顆粒平均半徑大小在長度方向與當前所建顆粒模型大小一致,標準差設為0.05,并通過Calculate Properties自動計算出顆粒質量、體積和轉動慣量X、Y、Z。

考慮到所用散體貯料為陶粒,顆粒數量多且在卸料過程中顆粒之間、顆粒與筒倉結構之間接觸時間較長,可通過加速度定律求解顆粒間的接觸碰撞力;同時顆粒之間、顆粒與筒倉結構之間無黏結,故選擇Hertz-Mindlin(no-slip)模型以模擬其接觸關系。依據材性試驗與相關文獻[13-14],顆粒接觸屬性參數與流動參數如表1所示。滿倉裝糧線為立筒倉頂部下降0.1 m,即貯料計算高度hn=1.18 m(hn為裝糧線至倉身下部的高度)。

表1 仿真參數Table 1 Simulation parameters

1.3 測點選定

根據立筒倉模型結構的對稱性和試驗條件,仿真模型數據測點的選取與試驗[12]所布測點位置保持一致,其中P1～P5用來采集貯料顆粒對倉壁產生的側壓力,分別對應測點深度0.2、0.5、0.8、1.1、1.256 m,且由于卸料方式為中心卸料,同一高度處的側壓力數據是相同的,故選取一列數據進行分析。因筒倉為三維對稱結構,故截取筒倉的中部垂直薄片進行數據分析,該垂直截面薄片設置為20 mm厚。在與側壓力測點P1～P5對應高度處,設置顆粒網格組(高度30 mm),其中V1～V5測點以獲取與側壓力測點高度一致的不同高度處顆粒速度,如圖2所示。

1.4 卸料工況

在進行卸料側壓力分析時,由于不同流速狀態下所產生的動態側壓力[15]與超壓系數會存在明顯不同,故將其分為兩種工況,流速慢與流速快工況分別命名為卸料a、卸料b,進而分析兩種流速卸料工況下糧食立筒倉靜態側壓力、動態側壓力與超壓系數。

2 模擬結果與分析

2.1 貯料流動狀態與速度分布

根據EDEM數值仿真和試驗結果,卸料流速不同時,兩種工況下的流動狀態性質相同,故進行流動狀態與速度分布分析時,選取卸料a(流速慢)工況為例,分別采用直接觀察和整體流量指數法進行卸料流動狀態分析。

2.1.1 直觀觀察卸料流態

通過截取不同時刻所選取的筒倉垂直截面薄片,直接觀察筒倉卸料流態,如圖3(a)～圖3(g)所示;同時根據筒倉卸料的實際工作情況,選取貯料頂部的5個顆粒進行標記,隱藏其它顆粒,觀察貯料顆粒的運動軌跡,如圖3(h)所示。

分圖題為卸料時間

由圖3(a)可知,當糧食達到穩定狀態后,筒倉開始進行卸料;由圖3(b)、圖3(c)可知,貯料為整體流動狀態,貯料上部表現為整體下降,下部呈漏斗流,卸料流動下降較緩慢;由圖3(d)、圖3(e)可知,當貯料下降至0.8 m測點深度處時,開始出現管狀流動趨勢,中間部位貯料流速明顯更快于倉壁兩側的貯料流速,倉壁附近的貯料顆粒向中心移動,并逐漸變為管狀流動階段,貯料下降明顯,流動均勻;由圖3(f)、圖3(g)可知,在卸料的末期,貯料表現出漏斗狀,流動加快,以卸料口正上方范圍內的顆粒速度最快先行流出,而靠近漏斗壁的顆粒由于阻礙作用形成死料區導致流出較晚,直至卸料完成。

隨著貯料顆粒的流動,在筒倉下部的顆粒之間、顆粒與倉壁之間的摩擦和擠壓更顯著,導致流動阻力增大,并在筒倉下部范圍內達到最大,但是此時筒倉上部的貯料顆粒不受此影響,阻力較小且流動較快,使得筒倉下部貯料更加密實,流動更加困難,相對應使得倉壁側壓力增加。

2.1.2 速度分布與整體流量指數

通過建立的V1～V5顆粒網格組,對整體卸料的貯料顆粒速度分布情況進行統計分析,根據模型內徑寬度,每個網格組在水平方向可均分為19個小網格,用于速度數據的采集,所提取速度代表為某一時刻顆粒經過該網格組時的豎向平均速度,如圖2(b)所示,在不同高度的豎向平均速度如表2所示。

表2 不同高度處豎向平均速度分布Table 2 Vertical average velocity distribution at different heights

通過整體流量指數[16](global flow index, MFI)對卸料流態情況進行量化分析,分析每個顆粒網格組速度分布,計算公式為

(1)

式(1)中:Vwall為筒倉內壁邊界處的貯料顆粒運動速度;Vcenter為筒倉內中心處的貯料顆粒運動速度。

當整體流量指數值大于0.3時該貯料流態為整體流,小于0.3時該貯料流態為漏斗流。因此根據表2計算不同高度處MFI值,其中V1和V2高度處的貯料表現為整體流,而從V3開始,貯料改變流態,逐漸呈漏斗流,結合直觀觀察卸料流態,并與試驗流態進行對比,兩者卸料時的流動狀態大致相同。

由圖4可知,在卸料整體過程中,中間貯料流動速度一直大于倉壁貯料流動速度,且筒倉倉壁[圖3(h)]豎向范圍內由上至下的不同深度位置的速度值逐漸減小,中間部位[圖3(h)]豎向范圍內由上至下的不同深度位置的速度值逐漸增大,即MFI值在逐漸減小,說明越靠近漏斗部分,貯料卸料流動狀態越接近漏斗流;在筒倉漏斗部分(V5測點處),因其接近于卸料口,故一直表現為漏斗流,筒倉在進行卸料作業時,其貯料所表現出的流態與所處于筒倉的位置有關。

括號中數值為MFI

2.2 立筒倉倉體側壓力分析

2.2.1 靜態側壓力

圖5為立筒倉滿倉靜態側壓力模擬值、試驗值和規范值對比結果。由圖5可知,靜態側壓力的模擬值同試驗值、規范值[17]一致,均隨著測點深度的增加均呈增大趨勢,最大值出現在倉壁與漏斗連接部位范圍內,靜態側壓力的模擬值較大于規范值,規范公式計算值具有一定的局限性且偏于保守,本次卸料仿真的模擬值與試驗值接近,為后續的動態側壓力與超壓系數分析提供可靠的依據。

圖5 立筒倉滿倉靜態側壓力模擬值、試驗值和規范值Fig.5 Simulation, test and specification values of static lateral pressure of silo full

2.2.2 動態側壓力與超壓系數

圖6為立筒倉兩種卸料工況下動態側壓力時程曲線。由圖6(a)可知,13 s時在0.8 m測點深度處壓力突增明顯,動態側壓力達8.79 kPa,是靜態側壓力的2.14倍;由圖6(b)可知,14 s時在0.8 m測點深度處壓力突增明顯,動態側壓力達8.77 kPa,是靜態側壓力的2.27倍。在卸料a與卸料b模擬時間0～150 s和0～48 s,貯料向下滑移,由于貯料顆粒之間、顆粒與倉壁之間的相互擠壓摩擦,一直對倉壁產生沖擊,致使倉壁上的P1～P4測點均出現了壓力突增的現象,并且發生震蕩現象;在漏斗處即測點P5處,由于離卸料口處最近,且在此處流態型式一直為漏斗流,會相應產生一部分死料區,故從卸料0 s開始,動態側壓力明顯迅速下降至某一段范圍內發生震蕩現象直至試驗結束,且無超壓現象發生。

圖6 立筒倉兩種卸料工況下動態側壓力時程曲線Fig.6 Dynamic lateral pressure time history curve of silo under two discharge conditions

在卸料過程中呈現出的震蕩現象,是由于卸料時貯料流動的過程中在1/3hn范圍內形成了超壓動力拱,使貯料內部反復的結拱和塌拱,拱腳對倉壁不斷進行沖擊造成。

由表3、圖7可知,EDEM仿真模擬所得的靜態、動態側壓力數據值與試驗做對比,兩者相差較小;筒倉倉壁且除漏斗處的各測點動態側壓力值相較于靜態側壓力值均有不同程度的增加,存在超壓現象,產生這種現象的原因是由于上部貯料重力和下部貯料的持續流出,導致發生成拱和破拱的現象,在0.8 m測點深度處成拱和破拱速度更快。在漏斗處即1.265 m測點深度發生動態側壓力減小的現象。隨著物料的流出,倉內貯料不斷地重復“結拱起始—結拱完成—拱塌落”的動態拱過程[18],因而倉壁卸料壓力呈現“較小—局部較大—較小”的震蕩分布。

表3 立筒倉側壓力模擬值與試驗值Table 3 Simulation and test values of silo lateral pressure

圖7 立筒倉靜態側壓力和動態側壓力Fig.7 Silo static lateral pressure and dynamic lateral pressure

由圖8可知,超壓系數整體呈現先升后降的趨勢,其中在0.8 m測點深度處最大,模擬卸料a超壓系數(2.14)和模擬卸料b超壓系數(2.27)均超過規范中水平側壓力修正系數Ch,其他測點超壓系數均在規范值內;在0.8 m測點深度處模擬卸料a、b相較于試驗卸料a、b的超壓系數誤差在3.4%、7.6%,差別較小且在合理誤差范圍內;且對比在不同測點深度處兩個工況超壓系數的大小,可發現不同卸料流速下,流速快工況的超壓系數均大于流速慢工況的超壓系數。

3 結論

通過對糧食立筒倉卸料進行EDEM(Engineering Discrete Element Method)仿真模擬分析,通過直接觀察和整體流量指數分析貯料流態,對比研究了不同流速下模擬值與試驗結果的貯料靜態壓力數據、動態側壓力分布規律與超壓系數,得出如下結論。

(1)筒倉下部1/3hn范圍內貯料的流動狀態發生顯著變化;貯料所表現出的流態與所處于筒倉的位置有關,越靠近漏斗部分,貯料卸料流動狀態越接近漏斗流;卸料的發展進程為整體流動—管狀流動—漏斗流。

(2)靜態側壓力數值模擬值與試驗值較為接近、動態側壓力分布規律一致,EDEM離散元卸料仿真模型動態側壓力研究具有可行性和合理性。

(3)卸料過程中,倉壁出現側壓力震蕩增大現象是由于在倉體下部1/3hn范圍內出現 “成拱—塌拱—再成拱—再塌拱”的動態拱過程。

(4)超壓系數最大值出現于0.8 m測點深度處;相同位置范圍內不同流速下的超壓系數均高于規范中水平側壓力修正系數Ch;不同測點深度處,流速快工況的超壓系數均大于流速慢工況的超壓系數。