基于EDEM的糧食立筒倉卸料動態(tài)側(cè)壓力研究

劉強(qiáng), 張茜茜, 劉濤, 丁永剛, 許啟鏗, 孫啟帥

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

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

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

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

1 離散元模型建立

1.1 立筒倉模型

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

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

1.2 貯料顆粒模型

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

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

表1 仿真參數(shù)Table 1 Simulation parameters

1.3 測點選定

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

1.4 卸料工況

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

2 模擬結(jié)果與分析

2.1 貯料流動狀態(tài)與速度分布

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

2.1.1 直觀觀察卸料流態(tài)

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

分圖題為卸料時間

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

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

2.1.2 速度分布與整體流量指數(shù)

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

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

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

(1)

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

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

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

括號中數(shù)值為MFI

2.2 立筒倉倉體側(cè)壓力分析

2.2.1 靜態(tài)側(cè)壓力

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

圖5 立筒倉滿倉靜態(tài)側(cè)壓力模擬值、試驗值和規(guī)范值Fig.5 Simulation, test and specification values of static lateral pressure of silo full

2.2.2 動態(tài)側(cè)壓力與超壓系數(shù)

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

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

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

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

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

圖7 立筒倉靜態(tài)側(cè)壓力和動態(tài)側(cè)壓力Fig.7 Silo static lateral pressure and dynamic lateral pressure

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

3 結(jié)論

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

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

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

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

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