基于離散元模擬筒倉貯料卸料成拱過程及筒倉壁壓力分布

劉克瑾，肖昭然，王世豪

?

基于離散元模擬筒倉貯料卸料成拱過程及筒倉壁壓力分布

劉克瑾，肖昭然，王世豪

（河南工業(yè)大學(xué)土木建筑學(xué)院，鄭州 450001）

筒倉卸料時貯料作用在倉壁上的卸料壓力出現(xiàn)驟然增大以及震蕩分布的現(xiàn)象，該文從貯料的散體顆粒性入手，采用離散元法和模型試驗法研究貯料在靜止儲糧狀態(tài)和卸料過程中的力學(xué)行為，從細(xì)觀顆粒層次探求卸料時貯料內(nèi)部土拱效應(yīng)與宏觀倉壁卸料壓力增大及產(chǎn)生震蕩的本質(zhì)聯(lián)系。模型為平底立筒倉，高1.0 m，寬0.5 m，卸料口直徑0.1 m，數(shù)值模型填充20 400個球形單元，模型試驗貯料為大豆。首先，通過分析卸料中倉底壓力分布的周期性變化規(guī)律，證實了卸料口附近拱效應(yīng)的存在。然后選取結(jié)拱起始、結(jié)拱完成及拱塌落3個時間點倉內(nèi)貯料的力鏈網(wǎng)絡(luò)、豎向應(yīng)力、橫向應(yīng)力、主應(yīng)力方向和速度場分布，分析了卸料時的拱效應(yīng)及其對倉壁卸料壓力分布的影響。研究發(fā)現(xiàn)，卸料中，筒倉底部的卸料口附近有拱形成，其跨度為卸料口直徑的4.0倍，高度為卸料口直徑的2.5倍。隨著物料的流出，卸料口附近的顆粒物質(zhì)遵循“拱形成-拱塌落”的動態(tài)規(guī)律，并據(jù)此提出了筒倉卸料的動態(tài)成拱機(jī)制。深高比0.35處，動態(tài)壓力修正系數(shù)最大為2.70。在深高比0.85處，結(jié)拱完成時的倉壁壓力達(dá)到峰值3.57 kPa。分析結(jié)果表明，拱的形成是倉壁壓力增大的根本原因，動態(tài)成拱機(jī)制則是宏觀倉壁壓力產(chǎn)生震蕩的根本原因，倉壁壓力峰值作用點和最大動態(tài)壓力修正系數(shù)作用點并不一致。動態(tài)成拱機(jī)制以及由此引發(fā)的倉壁卸料壓力分布規(guī)律，可為構(gòu)建機(jī)理研究的筒倉結(jié)構(gòu)安全設(shè)計提供參考。

應(yīng)力；數(shù)值方法；筒倉；倉壁壓力；拱；大豆；顆粒物質(zhì)；細(xì)觀研究

0 引 言

一定數(shù)量的顆粒聚集體具有復(fù)雜的力學(xué)特性，傳統(tǒng)的固體力學(xué)理論、流體力學(xué)理論和凝聚態(tài)物理很難對其做出準(zhǔn)確的解釋。筒倉作為貯存顆粒物質(zhì)的主要載體，其受力特性十分復(fù)雜。研究表明，筒倉散體物料的卸料過程中，倉壁受到的卸料壓力峰值可以達(dá)到靜態(tài)壓力的數(shù)倍[1-4]，大部分筒倉工程事故都是由這種過大的卸料壓力引起的。近年來，學(xué)者們對此進(jìn)行了探索并取得了初步的共識，即卸料過程中出現(xiàn)在卸料口附近的拱效應(yīng)是倉壁壓力增大的主要原因[5-7]。對于拱效應(yīng)的研究，由于傳統(tǒng)測量方法的局限性，室內(nèi)試驗僅能得到筒倉側(cè)壁及底部的邊界受力，對筒倉內(nèi)部散體顆粒體的受力特性不能進(jìn)行很好的觀測分析，因此也就無法直觀的觀察并量測到拱效應(yīng)。而近年來發(fā)展的先進(jìn)試驗方法如靜態(tài)光彈試驗（孫其誠等[8]）通過量測靜態(tài)堆積狀態(tài)下的顆粒間接觸力分布，發(fā)現(xiàn)了靜態(tài)堆積顆粒底部存在應(yīng)力凹陷現(xiàn)象，證明了拱效應(yīng)的存在。因卸料時顆粒物質(zhì)體系處于運動狀態(tài)，屬于大變形問題，由于等傾線的獲得需要進(jìn)行多次重復(fù)試驗，動態(tài)光彈應(yīng)力的定量分析及卸料的可重復(fù)性及精確性還有待提高[9]，不能對整個卸料過程進(jìn)行精確分析。傳統(tǒng)的室內(nèi)試驗無法解決的問題，從數(shù)值方法進(jìn)行探索是一條新思路。

當(dāng)前，數(shù)值方法中的有限元法和離散元法在模擬筒倉物料力學(xué)性能方面獲得了廣泛的應(yīng)用。尤其是離散元法，基于作用在顆粒體上的受力分析來確定單個顆粒體的運動行為，顆粒的運動遵循牛頓第二定律。與有限元法將顆粒集合體簡化為連續(xù)介質(zhì)不同，離散元法對問題的模擬是基于顆粒尺度（particle-level），可以從本質(zhì)上揭示顆粒集合體的宏觀力學(xué)行為。此外，離散元法還可以無損的觀察貯料內(nèi)部受力特性和流動規(guī)律，彌補模型試驗的不足。

學(xué)者們在用離散元法模擬筒倉的貯料和卸料過程，進(jìn)行了初步的探索，取得了有益的成果，然而也存在需要改進(jìn)的地方。這些研究多停留在宏觀的唯象學(xué)描述，并沒有充分發(fā)揮離散元法從細(xì)觀顆粒層次研究顆粒物質(zhì)力學(xué)行為的優(yōu)勢[10-13]。例如國內(nèi)學(xué)者們觀察到卸料中倉壁壓力出現(xiàn)“震蕩”變化，提出這一現(xiàn)象與“拱”效應(yīng)相關(guān)，然而其研究局限于現(xiàn)象的表述，并沒有深入的探討二者的本質(zhì)關(guān)系。再者，盡管有學(xué)者提出了拱效應(yīng)[14-17]，也推導(dǎo)了拱的形態(tài)與倉壁壓力的公式，然而這些推導(dǎo)均是基于連續(xù)介質(zhì)力學(xué)，并沒有從顆粒間的受力和傳力角度進(jìn)行分析，因而也沒有從本質(zhì)上揭示物料的宏觀力學(xué)表現(xiàn)。

鑒于此，本文采用離散元方法模擬筒倉貯料及中心卸料工況，借鑒顆粒物質(zhì)力學(xué)的基本理論，通過系統(tǒng)分析貯料的力鏈網(wǎng)絡(luò)分布、豎向應(yīng)力和橫向應(yīng)力分布、倉底和倉壁的卸料壓力分布，流速分布，主應(yīng)力方向分布，從細(xì)觀顆粒層次探討卸料中的拱效應(yīng)原理，包括拱的起始、形成和塌落，進(jìn)一步研究細(xì)觀的拱效應(yīng)如何影響宏觀的倉壁卸料壓力分布。最終從本質(zhì)上揭示貯料的宏觀力學(xué)行為。

1 模型建立

采用顆粒離散元軟件PFC模擬貯料及中心卸料工況，并制作了相同尺寸筒倉模型的室內(nèi)貯料及卸料試驗。模型筒倉高1 m，直徑0.5 m，卸料口直徑0.1 m，位于筒倉底部中心，模型縮放比為1:20。倉體底部測墻長度為0.05 m，從左到右編號依次為1～8。倉體左右兩側(cè)各設(shè)置倉壁壓力監(jiān)測墻10個，測墻高度均為0.1m，其中左側(cè)測墻沿高度方向從下到上編號依次為11～19；右側(cè)測墻沿高度方向從下到上編號依次為21～29。室內(nèi)模型試驗中，筒倉形狀為圓柱形，內(nèi)徑0.5 m，高1 m，沿倉壁兩側(cè)母線對稱布置壓力計，每側(cè)6個共計12個（圖1b），測量倉壁壓力。試驗采用DYB-3型電阻應(yīng)變式壓力計，其量程為20 kPa，分辨率≤0.083% F.S，綜合誤差＜0.8% F.S，數(shù)據(jù)采集儀采用DH5922動態(tài)應(yīng)變分析測試系統(tǒng)。共進(jìn)行了5次平行試驗，每次試驗取左右兩側(cè)壓力計均值，文中給出的貯料壓力及卸料壓力為5次試驗值的均值。模型筒倉及室內(nèi)試驗?zāi)Ｐ腿鐖D1所示。

注：1~8，11~30表示測墻編號，下同。

本文模擬的貯料為豆類等似球形散粒體，采用圓球單元，參考前人對大豆籽粒的彈性模量研究成果[18]，設(shè)置顆粒的剛度參數(shù)（表1）。依據(jù)大豆與倉壁材料摩擦系數(shù)研究成果[19]，設(shè)定墻摩擦系數(shù)為0.4，該系數(shù)可以代表墻體材料為鋼板或混凝土。大豆的內(nèi)摩擦角依據(jù)大豆內(nèi)摩擦角研究成果[20]，綜合考慮各項影響因素，選為32°，對應(yīng)顆粒的摩擦系數(shù)為0.6。大豆顆粒密度與頂部壓力有關(guān)，本文中，經(jīng)Janssen公式計算大豆的頂部壓力分布為0～1.80 kPa，參考大豆堆壓縮密度與體變模量研究成果[21]設(shè)置球形顆粒單元的顆粒密度為718 kg/m3，具體顆粒單元參數(shù)示于表1。

表1 筒倉離散元法模型材料參數(shù)

2 貯料工況分析

數(shù)值模型中，倉內(nèi)貯料采用點源法生成，當(dāng)顆粒全部下落完成后，測墻輸出倉壁及倉底的壓力曲線呈水平直線時認(rèn)為模型達(dá)到平衡狀態(tài)，定義為貯料工況，此時模型倉中顆粒共計20 400個。室內(nèi)試驗?zāi)Ｐ椭?，滿倉后靜置1 h，土壓力盒讀數(shù)穩(wěn)定，定義為貯料工況。貯料工況下的倉壁壓力分布與Janssen公式[22]計算值對比如圖2所示。本文中，深高比為貯料深度與筒倉高度之比。

圖2 倉壁側(cè)壓力模擬值、Jassen解和試驗值

圖2顯示數(shù)值模擬結(jié)果、試驗結(jié)果與工程界普遍認(rèn)可的Janssen公式計算值有較好的一致性。以試驗值為參考值，模擬值和Jassen值與試驗值的最大誤差不超0.04和0.09，證明了離散元法的可靠性。為了觀察貯料工況下顆粒接觸力的分布情況，用PFC軟件輸出貯料工況下倉內(nèi)顆粒物質(zhì)的力鏈分布、孔隙率分布、豎向應(yīng)力分布和橫向應(yīng)力分布，如圖3所示。

圖3a中，顆粒間接觸力由線段表示，線段越粗，表示顆粒間接觸力越大，力鏈強度越大。由圖3a可以看出，貯料工況下，筒倉下部力鏈強度大且分布較為密集，相互鏈接形成架拱[23]。筒倉上部力鏈強度小且分布較為稀疏。筒倉底部，力鏈形成由實線標(biāo)注的以及其他大小不一的弧形拱狀結(jié)構(gòu)，這與房營光等[24]發(fā)現(xiàn)的試驗現(xiàn)象一致，倉底接觸部分力鏈在水平方向上分布的不均勻性是壓力凹陷[8]的主要原因。由圖3b可以看出，由于填料時物料的隨機(jī)排布特性，筒倉內(nèi)的孔隙率分布并不均勻，存在類似虛線標(biāo)注區(qū)域等孔隙率極小的顆粒塊體。圖3c中，豎向應(yīng)力在豎直方向上由上到下逐漸增大，在同一高度，豎向應(yīng)力在水平方向上的分布并不均勻，這種分布在筒倉中上部尤為明顯。圖3d中，橫向應(yīng)力沿豎直方向由上到下非線性增大。

圖3 貯料靜工況顆粒物質(zhì)的力鏈網(wǎng)絡(luò)、孔隙率、豎向應(yīng)力和橫向應(yīng)力分布

對比圖3a和3b發(fā)現(xiàn)，在顆粒塊體內(nèi)部，力鏈強度較弱且分布密集。弱力鏈在整個筒倉內(nèi)都有分布，而強力鏈主要分布在孔隙率較小的區(qū)域。對比圖3a、3b和3c發(fā)現(xiàn)，在豎向應(yīng)力較大的區(qū)域內(nèi)，強度較大的力鏈方向與重力方向夾角較?。辉跈M向應(yīng)力較大的區(qū)域內(nèi)，強度較大的力鏈方向與重力方向夾角較大。說明強度較大的力鏈的方向影響區(qū)域內(nèi)豎向應(yīng)力及橫向應(yīng)力的分布。

3 卸料過程

在分析卸料過程之前，參考已有研究結(jié)論[13-14,16]：貯料卸料流動時的結(jié)拱－破拱引發(fā)了筒倉卸料時倉壁側(cè)壓力不斷的上下波動，拱效應(yīng)位于卸料口附近。本文嘗試觀察倉底壓力的分布，預(yù)判是否存在拱效應(yīng)。若不存在拱效應(yīng)，則卸料中倉底壓力分布可能是單調(diào)增大、不變或減?。蝗舸嬖诠靶?yīng)，即卸料口附近有拱的形成和塌落，那么拱的形成會使位于拱腳部分的倉底壓力增大，拱腳兩側(cè)的倉底壓力則會減小。

3.1 倉底壓力分析

自卸料開始，記錄各個測墻的壓力變化，倉底測墻編號及位置見圖4，其中1、2、7、8號測墻遠(yuǎn)離卸料口，3、4、5、6號測墻靠近卸料口。

圖4 拱腳位置

研究拱效應(yīng)時，考慮到顆粒物質(zhì)生成的隨機(jī)性，歷次結(jié)拱的拱腳位置可能會有差別，因此，需進(jìn)行多次模擬取統(tǒng)計學(xué)分析結(jié)果。本文計算時做了5次數(shù)值模擬。倉底壓力的谷值和和峰值變化趨勢以及倉底的壓力分布基本一致，誤差不超過均值的5%。分析后發(fā)現(xiàn)結(jié)拱的拱腳位置均處于1、2號和7、8號測墻，限于篇幅此處選取有代表性的倉底壓力分布圖5進(jìn)行分析。

圖5為卸料初期，且卸料口附近顆粒速度穩(wěn)定時，作用在倉底上的卸料壓力分布圖。由圖5可知：

1）卸料過程中倉底及倉壁的壓力一直處于震蕩變化狀態(tài)，這一現(xiàn)象與文獻(xiàn)[13,25-26]的研究相符?？拷读峡诘膫}底壓力較小，5號測墻倉底壓力與倉底峰值壓力之比不超過0.41，4號測墻不超過0.39。遠(yuǎn)離卸料口倉底壓力較大，如1、2號和7、8號測墻，倉底壓力與倉底峰值壓力之比在0.5與1.0之間周期性變化。

2）倉底的壓力變化具有一定的周期性。這一現(xiàn)象證明了卸料口附近拱效應(yīng)的存在。尤其是1、2號和7、8號周期性壓力變化明顯，且其承擔(dān)的壓力值較大。根據(jù)拱的存在會使拱腳部分的倉底壓力增大這一物理現(xiàn)象，進(jìn)行了5次數(shù)值模擬，進(jìn)行統(tǒng)計分析后發(fā)現(xiàn)，1、2號和7、8號是拱腳的位置。

3）研究拱效應(yīng)的節(jié)點選擇。研究表明，拱的形成會使局部（拱腳部分）倉底壓力增大，拱的塌落會使局部倉底壓力減小[27-28]。以此為依據(jù)：在壓力始終較大且隨時間變化明顯的1、2號和7、8號底墻的壓力變化曲線上，選取極小值點A作為結(jié)拱的起始點，選取極大值點B作為結(jié)拱的完成點，選取B點之后迅速減小后的一點C作為拱塌落的代表點。在此基礎(chǔ)上，分別研究對應(yīng)點A、B和C時刻的筒倉內(nèi)部力鏈網(wǎng)絡(luò)分布，豎向應(yīng)力分布，橫向應(yīng)力分布及顆粒速度分布，并據(jù)此從細(xì)觀顆粒層次分析拱的形成－塌落機(jī)理及其對倉內(nèi)貯料應(yīng)力分布的影響。

注：A、B、C分別為結(jié)拱起始點、結(jié)拱完成點和拱塌落點

3.2 細(xì)觀顆粒層次的拱效應(yīng)分析

3.2.1 結(jié)拱起始

圖6a為結(jié)拱起始時刻點A的貯料力鏈圖，圖6b為豎向應(yīng)力分布，圖6c為橫向應(yīng)力（也即倉壁壓力）分布，圖6d為顆粒速度。

根據(jù)力鏈稀疏程度、力的傳力角度、力鏈強弱、流速不同，可以將倉內(nèi)顆粒體系劃分為4個區(qū)域見圖6a中區(qū)域1，區(qū)域2，區(qū)域3，區(qū)域4，觀察圖6可知：

1）區(qū)域1和4內(nèi)的力鏈較為稀疏，且區(qū)域1內(nèi)的長度較大的力鏈方向與豎向夾角較大；強力鏈主要分布在區(qū)域2和3內(nèi)。相應(yīng)的區(qū)域1和4內(nèi)豎向應(yīng)力和橫向應(yīng)力均較小，區(qū)域2和3內(nèi)有較大的豎向應(yīng)力和橫向應(yīng)力分布。

2）圖6d顯示，區(qū)域4內(nèi)的顆粒流速較其他區(qū)域更大，區(qū)域3頂部以上的顆粒流動直徑約為筒倉直徑的3/4，在區(qū)域2和3底部附近，越靠近倉底，流動半徑越小，即貯料的流態(tài)為漏斗流。

分析可知：在強度較大的力鏈分布的區(qū)域內(nèi)，豎向應(yīng)力或橫向應(yīng)力較大，這是因為應(yīng)力的分布是由區(qū)域內(nèi)力鏈強度決定的。對比圖6b、6c和6d發(fā)現(xiàn)，橫向應(yīng)力的分布與區(qū)域內(nèi)顆粒速度無明顯聯(lián)系；而在豎向應(yīng)力較大的區(qū)域內(nèi)，顆粒速度較小，反之亦然，說明區(qū)域內(nèi)顆粒速度主要由豎向應(yīng)力的分布決定。

注：區(qū)域1、2、3、4分別表示力鏈稀疏區(qū)、力鏈密集區(qū)、成拱區(qū)和快速流動區(qū)。

3.2.2 結(jié)拱完成

圖7為結(jié)拱完成時刻點B的貯料力鏈圖，豎向應(yīng)力分布圖，橫向應(yīng)力分布圖及顆粒速度圖。

觀察圖7可以發(fā)現(xiàn)：

1）區(qū)域1和2中均有強力鏈分布，倉壁附近力鏈以豎向為主，區(qū)域3中有力鏈呈拱形分布，區(qū)域4中力鏈強度較弱，且分布稀疏。相應(yīng)的，區(qū)域1和2中均有較大的豎向應(yīng)力和橫向應(yīng)力分布，較大的豎向應(yīng)力主要集中分布在倉壁及區(qū)域2底部附近；區(qū)域3內(nèi)有較大的豎向應(yīng)力和橫向應(yīng)力呈拱形分布，符合拱的受力特性。

2）圖7d顯示，區(qū)域1內(nèi)顆粒速度很小，頂部的顆粒速度為0，停止了流動，正是由于拱的形成使得拱的形成區(qū)域3及其上方流動區(qū)域的顆粒運動速度減小。

對比圖7a、7b、7c和7d發(fā)現(xiàn)，區(qū)域1內(nèi)強度較大的力鏈主要影響橫向應(yīng)力的大小及其分布，區(qū)域2內(nèi)強度較大的力鏈主要影響豎向應(yīng)力的大小及其分布，說明區(qū)域1內(nèi)強度較大的力鏈方向以橫向為主，而區(qū)域2內(nèi)強度較大的力鏈方向以豎向為主。區(qū)域3內(nèi)形成了跨度約為卸料口直徑4倍，高度約為卸料口直徑2.5倍的拱結(jié)構(gòu)。區(qū)域4位于拱圈以下，該區(qū)遍布弱力鏈，豎向應(yīng)力及橫向應(yīng)力接近零，相應(yīng)的該區(qū)顆粒的流動速度也最大。

圖7 結(jié)拱完成時顆粒物質(zhì)的力鏈網(wǎng)絡(luò)、豎向應(yīng)力、橫向應(yīng)力和速度分布

從圖6結(jié)拱開始到圖7結(jié)拱完成，拱的形成主要來自于區(qū)域3內(nèi)力鏈強度由小到大的變化。由于區(qū)域3內(nèi)拱的形成，阻礙了區(qū)域2和區(qū)域1內(nèi)顆粒的流動（直觀表現(xiàn)是區(qū)域1和區(qū)域2內(nèi)顆粒速度大幅度減?。瑓^(qū)域1和2內(nèi)的豎向應(yīng)力和橫向應(yīng)力均有明顯增加，且橫向應(yīng)力增加趨勢更為明顯，此時，橫向應(yīng)力向兩側(cè)傳遞作用在筒倉側(cè)壁的橫向壓應(yīng)力必然也會增大，由此證明了拱的形成導(dǎo)致卸料壓力增大。該結(jié)論與已發(fā)表的土力學(xué)中的拱效應(yīng)研究成果[28]相容。

3.2.3 拱塌落

圖8為拱塌落時刻點C的貯料力鏈網(wǎng)絡(luò)分布圖，豎向應(yīng)力分布圖，橫向應(yīng)力分布圖及顆粒速度圖。

圖8a和8b顯示，強度較大的力鏈及豎向應(yīng)力主要分布于區(qū)域1中下部、區(qū)域2倉壁附近以及區(qū)域3底部。圖8c顯示，區(qū)域內(nèi)中下部橫向應(yīng)力較大，區(qū)域2內(nèi)橫向應(yīng)力次之，而區(qū)域3和4內(nèi)橫向應(yīng)力接近于0。圖8d顯示，區(qū)域4內(nèi)的顆粒速度較其他區(qū)域更大。

由圖7結(jié)拱完成到圖8拱塌落，貯料顆粒體系的力鏈強度由大變小，區(qū)域2和3內(nèi)的力鏈強度減小尤為明顯。區(qū)域1、2和3內(nèi)的豎向應(yīng)力和橫向應(yīng)力也有較大程度的減小，說明拱處于塌落狀態(tài)。以此可以得出，拱塌落時刻點C的倉壁壓力小于結(jié)拱完成時刻點B的倉壁壓力。

圖8 拱塌落時顆粒物質(zhì)的力鏈網(wǎng)絡(luò)、豎向應(yīng)力、橫向應(yīng)力和速度分布

3.3 主應(yīng)力方向分析

PFC采用Christoffersen[29]和Fortin[30]提出的顆粒間接觸力計算范圍內(nèi)應(yīng)力的方法計算一個圓形區(qū)域內(nèi)的應(yīng)力張量，為更精確的判斷拱的形成對倉內(nèi)應(yīng)力分布的影響，并借鑒賴漢江等[31]研究成果，本文將筒倉劃分為200個大小相同的正方形區(qū)域，每個區(qū)域內(nèi)均布置有內(nèi)切于正方形的圓形測量圓，橫向應(yīng)力σ、切向應(yīng)力σ及豎向應(yīng)力σ由測量圓測定，并通過式（3）求得主應(yīng)力及絕對值最大的主應(yīng)力方向。將貯料工況及結(jié)拱起始時刻點A，結(jié)拱完成時刻點B和拱塌落時刻點C的筒倉內(nèi)絕對值最大主應(yīng)力方向（以下簡稱主應(yīng)力方向）分布示于圖9。

圖9 不同時刻主應(yīng)力方向分布圖

圖9a顯示，貯料工況下主應(yīng)力方向以豎向為主，且除筒倉頂部主應(yīng)力方向較雜亂外，其他區(qū)域內(nèi)主應(yīng)力方向分布有較好的規(guī)律性：距倉壁越近，主應(yīng)力方向與重力方向夾角越大。圖9b顯示，結(jié)拱起始時刻，筒倉下部側(cè)壁附近主應(yīng)力方向以豎向為主，但其他區(qū)域內(nèi)的主應(yīng)力方向較貯料工況發(fā)生了較大偏轉(zhuǎn)。圖9c顯示，結(jié)拱完成時刻，在懸鏈線處形成了弧形拱狀結(jié)構(gòu)，拱下方的主應(yīng)力方向分布較為雜亂，在拱頂?shù)恼戏?，由于拱的形成阻礙了上方顆粒的流動，對上方顆粒起到了支撐作用，拱頂上方的主應(yīng)力方向以豎向為主。圖9d拱塌落時刻，除筒倉下部側(cè)壁附近主應(yīng)力方向以豎向為主，有著較好的規(guī)律性外，其他區(qū)域內(nèi)主應(yīng)力方向分布較為雜亂。

對比圖9b、9c和9d發(fā)現(xiàn)，雖然圖9b和9d中懸鏈線內(nèi)部主應(yīng)力方向有拱形分布，但其附近的主應(yīng)力方向分布較為雜亂，說明懸鏈線附近沒有拱的形成，而圖9c中懸鏈線附近的主應(yīng)力方向分布較圖9b和9d有更好的規(guī)律性，說明懸鏈線附近有拱的形成。在拱的形成過程中，由于拱對上方顆粒的支撐作用，拱頂上方的主應(yīng)力方向與重力方向的夾角減?。辉诠暗乃溥^程中，由于失去了拱的支撐作用，拱頂上方的主應(yīng)力方向發(fā)生偏轉(zhuǎn)，與重力方向的夾角增大。

3.4 宏觀倉壁卸料壓力分析

為明確拱的形成與塌落對筒倉立壁（下稱倉壁）壓力分布的影響，繪制貯料狀態(tài)及結(jié)拱起始時刻點A、結(jié)拱完成時刻點B及拱塌落時刻點C的倉壁壓力分布如圖10a。值得一提的是，由于是中心卸料，原則上講左右測墻壓力應(yīng)該是對稱分布，然而由于離散元法在生成顆粒單元時，球體單元排布的隨機(jī)性，以及其散體顆粒性質(zhì)，左右測墻的壓力分布并非完全對稱，而是有微小的差別，以左測墻為參考，這種差別不超過8%。圖10a中倉壁側(cè)壓力采用的是左右測墻壓力的平均值。

參考國家規(guī)范《糧食鋼板筒倉設(shè)計規(guī)范》（GB50322-2011）第4.2.3條，定義動態(tài)壓力修正系數(shù)（overpressure coefficient）為卸料中動態(tài)倉壁壓力與靜態(tài)倉壁壓力之比，圖10b給出了動態(tài)壓力修正系數(shù)隨深高比的分布圖。

圖10 倉壁側(cè)壓力分布

觀察圖10，在貯料工況、結(jié)拱起始點、結(jié)拱完成點及拱塌落時刻，倉壁壓力均隨著深度的增加非線性增大。結(jié)拱起始時刻的倉壁壓力較貯料靜壓力相差不大。在筒倉的淺層（深高比小于0.2），動態(tài)壓力修正系數(shù)小于1，即動態(tài)倉壁壓力小于靜態(tài)倉壁壓力，分析這部分的工程實際意義不大。下面著重分析工程界關(guān)心的深層倉壁壓力分布。

為了驗證模擬結(jié)果，圖10a中附上按照國家規(guī)范[32]計算的倉壁卸料壓力值與本文室內(nèi)試驗結(jié)果及Hala中心卸料的模型試驗研究成果[33]。在筒倉中下部，也即深高比大于0.45，結(jié)拱完成值超過了規(guī)范計算值，用規(guī)范預(yù)測本文的倉壁壓力值可能是危險的。結(jié)拱完成值要低于Hala的試驗值，這可能是由于Hala在試驗中采用貯料為玉米，本文采用的為圓球顆粒，圓球顆粒的粒間咬合力很小，這一現(xiàn)象與Koby?ka[34]在做數(shù)值模擬與驗證時的發(fā)現(xiàn)相似。從整體上看Hala的室內(nèi)模型試驗結(jié)果、本文的試驗結(jié)果與模擬值（結(jié)拱完成值）有較好的一致性，即二者均存在卸料壓力的震蕩現(xiàn)象，且卸料壓力峰值均出現(xiàn)在筒倉下部深高比0.8～0.85之間。

深高比超過0.2時，結(jié)拱完成時的超壓系數(shù)大于1。深高比0.35處，動態(tài)壓力修正系數(shù)最大為2.70（圖10b中的b點）。在深高比0.85處，結(jié)拱完成時的倉壁壓力最大，達(dá)到3.57 kPa。這一結(jié)論與本文室內(nèi)試驗結(jié)果及Hala[33]的室內(nèi)模型試驗結(jié)果接近，此處的超壓系數(shù)為2.18（圖10b中的a點）。即倉壁壓力峰值作用點和最大動態(tài)壓力修正系數(shù)作用點并不一致。

3.5 細(xì)觀動態(tài)成拱機(jī)制對宏觀倉壁卸料壓力分布的影響

綜合離散元分析結(jié)果，結(jié)合圖6－圖10可知，卸料中的成拱機(jī)制是動態(tài)的，拱效應(yīng)的存在是卸料中倉壁壓力增大以及震蕩變化的主要原因。

卸料是顆粒流出卸料口的過程，在拱形成之前，位于倉體上部的顆粒體系遍布弱力鏈，其顆粒流速較大，宏觀表現(xiàn)是此時倉壁卸料壓力較小。由于顆粒流動時重新排列組合，其間必然存在相對位移產(chǎn)生摩擦阻力，這些阻力沿著流動方向逐漸增大，在卸料口處達(dá)到最大值。此時，筒倉上部阻力小，顆粒流速較快（圖6d），而下部流速較慢，上部顆粒的流動會對下部顆粒產(chǎn)生壓實作用，從而導(dǎo)致卸料口附近的顆粒物質(zhì)擠壓成拱。成拱的細(xì)觀表現(xiàn)是料口附近力鏈強度增大、豎向應(yīng)力及橫向應(yīng)力的增大，主應(yīng)力方向的規(guī)律化，宏觀表現(xiàn)是拱腳處的倉底壓力達(dá)到局部極大值，作用在倉壁上的卸料壓力增大，最大達(dá)到靜態(tài)貯料倉壁壓力的2.7倍。

由于拱的形成，阻礙了上部顆粒的運動，細(xì)觀表現(xiàn)是拱頂周圍顆粒流速減小，力鏈強度增大，豎向應(yīng)力和橫向應(yīng)力的增大，拱圈內(nèi)部的顆粒流速則達(dá)到最大，遍布弱力鏈。一方面，拱圈內(nèi)部顆粒不斷的流出，創(chuàng)造著新的臨空面。另一方面，拱頂上方的顆粒由于受到拱的阻礙，施加給拱更大的擠壓力，當(dāng)擠壓力累加到拱無法承受時，拱開始塌落，此時，拱體的力鏈強度迅速減小，豎向應(yīng)力和橫向應(yīng)力減小，主應(yīng)力方向的無序化，宏觀表現(xiàn)則是倉壁卸料壓力的減小。

之后，隨著物料的流出，倉內(nèi)貯料不斷的重復(fù)“結(jié)拱起始-結(jié)拱完成-拱塌落”的動態(tài)成拱過程，因而反映在宏觀上則是倉壁卸料壓力呈現(xiàn)“較小-局部極大-較小”的震蕩分布。

4 結(jié) 論

本文采用離散元方法和室內(nèi)模型試驗?zāi)M了筒倉貯料工況及卸料工況，模型試驗結(jié)果驗證了離散元法結(jié)果。通過系統(tǒng)的研究結(jié)拱起始、結(jié)拱完成和拱塌落等時間節(jié)點的貯料力鏈網(wǎng)絡(luò)、豎向應(yīng)力和橫向應(yīng)力、顆粒速度及倉壁壓力分布，從細(xì)觀顆粒層次探討了卸料過程中的動態(tài)成拱機(jī)制，以及由此引發(fā)的宏觀倉壁側(cè)壓力分布。主要得出以下結(jié)論：

1）貯料工況下，同一高度的豎向應(yīng)力在水平方向上非線性變化，同一橫向位置的橫向應(yīng)力在豎直方向上非線性變化。

2）卸料中，倉壁及倉底壓力處于周期性變化的狀態(tài)。

3）卸料中，筒倉底部的卸料口附近有拱形成，其跨度約為卸料口直徑的4倍，高度約為卸料口直徑的2.5倍。

4）拱的形成導(dǎo)致卸料過程中倉壁壓力增大，拱塌落使得倉壁壓力減小。動態(tài)成拱機(jī)制即“結(jié)拱（結(jié)拱起始-結(jié)拱完成）-拱塌落”導(dǎo)致倉壁的卸料壓力呈現(xiàn)“較小-局部極大-較小”的震蕩分布。

5）工程界關(guān)心的倉壁壓力峰值作用點在筒倉下部深高比0.85處。

6）本文的動態(tài)壓力修正系數(shù)最大值為2.70，位于深高比0.35處。

[1] 陳長冰，梁醒培. 筒倉卸料過程的離散元模擬分析[J]. 糧油食品科技，2008，16(1)：11－13.

Chen Changbing, Liang Xingpei. Analysis on discrete element simulation of discharging in silos[J]. Science and Technology of Cereals, Oils and Foods, 2008, 16(1): 11－13. (in Chinese with English abstract)

[2] Koby?ka R, Molenda M. DEM simulations of loads on obstruction attached to the wall of a model grain silo and of flow disturbance around the obstruction[J]. Powder Technology, 2014, 256: 210－216.

[3] 劉定華，郭明昌，楊建斌. 筒倉倉壁動態(tài)壓力的測試和分析[J]. 西安建筑科技大學(xué)學(xué)報：自然科學(xué)版，1989，21(1)：108－113.

Liu Dinghua, Guo Mingchang, Yang Jianbin. Test and analysis on dynamical pressure to the wall of silo[J]. Journal of Xi’an University of Architecture & Technology: Natural Science Edition, 1989, 21(1): 108－113. (in Chinese with English abstract)

[4] 李國柱，劉定華. 筒倉動態(tài)壓力的計算和測試[J]. 寧波高等?？茖W(xué)校學(xué)報，2000，12(4)：1－7.

Li Guozhu, Liu Dinghua. Calculation of dynamic pressure on reinforced concrete silo[J]. Journal of Ningbo College, 2000, 12(4): 1－7. (in Chinese with English abstract)

[5] Wang Peijun, Zhu Xulin, Liu Mei, et al. Buckling behaviors and simplified design method for steel silos under locally distributed axial load[J]. Journal of Constructional Steel Research, 2017, 134: 114－134.

[6] 趙松. 筒倉貯料壓力分析及其應(yīng)用[D]. 武漢：武漢理工大學(xué)，2013.

Zhao Song. Analysis and Application for Silo Storage Pressure[D]. Wuhan: Wuhan University of Technology, 2013. (in Chinese with English abstract)

[7] 車穎文. 筒倉的受力情況和穩(wěn)定性研究[D]. 武漢：武漢理工大學(xué)，2011.

Che Yingwen. Research on the Stress Behavior and Stability of the Silos[D]. Wuhan: Wuhan University of Technology, 2013. (in Chinese with English abstract)

[8] 孫其誠，王光謙. 靜態(tài)堆積顆粒中的力鏈分布[J]. 物理學(xué)報，2008，57(8)：4667－4674.

Sun Qicheng, Wang Guangqian. Force distribution in static granular matter in two dimensions[J]. Acta Physica Sinica, 2008, 57(8): 4667－4674. (in Chinese with English abstract)

[9] 鄒同彬，陸渝生，連志穎. 動光彈實驗中的主應(yīng)力分離方法[J]. 解放軍理工大學(xué)學(xué)報：自然科學(xué)版，2004，5(4)：72－76.

Zou Tongbin, Lu Yusheng, Lian Zhiying. Separation of principal stress in dynamic photoelastic experiment[J]. Journal of PLA University of Science and Technology, 2004, 5(4): 72－76. (in Chinese with English abstract)

[10] 譚援強，肖湘武，鄭軍輝，等. 錐形改流體下部孔徑對筒倉卸料流態(tài)的影響[J]. 農(nóng)業(yè)工程學(xué)報，2016，32(19)：82－87.

Tan Yuanqiang, Xiao Xiangwu, Zheng Junhui, et al. Effect of outlet diameter of cone-in-cone insert on silo flow pattern[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(19): 82－87. (in Chinese with English abstract)

[11] 俞良群，邢紀(jì)波. 筒倉裝卸料時力場及流場的離散單元法模擬[J]. 農(nóng)業(yè)工程學(xué)報，2000，16(4)：15－19.

Yu Liangqun, Xing Jibo. Discrete element method simulation of forces and flow fields during filling and discharging materials in silos[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2000, 16(4): 15－19. (in Chinese with English abstract)

[12] 肖昭然，王軍，何迎春. 筒倉側(cè)壓力的離散元數(shù)值模擬[J]. 河南工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2006，27(2)：10－12，16.

Xiao Zhaoran, Wang Jun, He Yingchun. The DEM simulation of the lateral pressure of silo[J]. Journal of Henan University of Technology: Natural Science Edition, 2006, 27(2): 10－12, 16. (in Chinese with English abstract)

[13] 張昭，劉克瑾，肖昭然，等. 筒倉倉壁摩擦對倉壁側(cè)壓力影響的研究[J]. 河南工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2017，38(5)：88－92.

Zhang Zhao, Liu Kejin, Xiao Zhaoran, et al. Influence of the wall friction on the lateral pressure of silo wall[J]. Journal of Henan University of Technology: Natural Science Edition, 2017, 38(5): 88－92. (in Chinese with English abstract)

[14] Oldal I, Keppler I, Csizmadia B, et al. Outflow properties of silos: The effect of arching[J]. Advanced Powder Technology, 2012, 23(3): 290－297.

[15] Goda T J, Ebert F. Three-dimensional discrete element simulations in hoppers and silos[J]. Powder Technology, 2005, 158(1/2/3): 58－68.

[16] Matchett A J. The shape of the cohesive arch in hoppers and silos-Some theoretical considerations[J]. Powder Technology, 2007, 171(3): 133－145.

[17] Cheng Y M, Chau K T, Xiao L J, et al. Flow pattern for a silo with two layers of materials with single or double openings[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(9): 1278－1286.

[18] 唐福元，馮家暢，程緒鐸. 大豆籽粒彈性模量的測定與研究[J]. 糧食儲藏，2015，44(5)：6－9.

Tan Fuyuan, Feng Jiachang, Cheng Xuduo. Determination and study on the elastic modulus of soybean seeds[J]. Grain Storage, 2015, 44(5): 6－9. (in Chinese with English abstract)

[19] 馮家暢，程緒鐸，嚴(yán)曉婕. 大豆與倉壁材料摩擦系數(shù)的研究[J]. 大豆科學(xué)，2014，33(5)：787－789+794.

Feng Jiachang, Cheng Xuduo, Yan Xiaojie. Research on the coefficient of friction between soybean and silo wall material[J]. Soybean Science, 2014, 33(5): 787－789+794. (in Chinese with English abstract)

[20] 程緒鐸，陸琳琳，石翠霞，等. 大豆內(nèi)摩擦角的測定與實驗研究[J]. 糧食儲藏，2010，39(5)：12－15.

Cheng Xuduo, Lu Linlin, Shi Cuixia, et al. Measurement and experimenta on internal friction angle of soybean[J]. Grain Storage, 2010, 39(5): 12－15. (in Chinese with English abstract)

[21] 馮家暢，程緒鐸，杜小翠，等. 大豆堆壓縮密度與體變模量研究[J]. 中國糧油學(xué)報，2016，31(12)：112－117.

Feng Jiachang, Cheng Xuduo, Du Xiaocui, et al. Compressive density and bulk modulus of soybean[J]. Journal of the Chinese Cereaks and Oils Association, 2016, 31(12): 112－117. (in Chinese with English abstract)

[22] Janssen H A. Experiments about pressure of grain in silos[J]. VDI, 1895: 1045－1049.

[23] 孫其誠，辛海麗，劉建國，等. 顆粒體系中的骨架及力鏈網(wǎng)絡(luò)[J]. 巖土力學(xué)，2009，30(增刊1)：83－87.

Sun Qicheng, Xin Haili, Liu Jianguo, et al. Skeleton and force chain network in static granular material[J]. Rock and Soil Mechanics, 2009, 30(Supp.1): 83－87. (in Chinese with English abstract)

[24] 房營光，侯明勛，谷任國，等. 樁承式路堤中土拱效應(yīng)產(chǎn)生過程可視化分析[J]. 巖土工程學(xué)報，2015，37(9)：1678－1684.

Fang Yingguang, Hou Mingxun, Gu Renguo, et al. Visual analysis of initiation of soil arching effect in piled embankments[J]. Rock and Soil Mechanics, 2015, 37(9): 1678－1684. (in Chinese with English abstract)

[25] 程奇鵬，孫巍巍，盧賽. 基于PFC3D的淺圓倉偏心卸料離散元研究[J]. 土木工程與管理學(xué)報，2016，33(6)：43－47.

Cheng Qipeng, Sun Weiwei, Lu Sai. Discrete element analysis of squat silo under eccentric discharge by PFC3D[J]. Journal of Civil Engineering and Management, 2016, 33(6): 43－47. (in Chinese with English abstract)

[26] 張大英，許啟鏗，王樹明，等. 筒倉動態(tài)卸料過程側(cè)壓力模擬與驗證[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(5)：272－278，316.

Zhang Daying, Xu Qikeng, Wang Shuming, et al. Simulation and experimental validation of silo wall pressure during discharging[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(5): 272－278, 316. (in Chinese with English abstract)

[27] Cannavaccuiolo A, Barletta D, Dosi G, et al. Arch-Free flow in aerated silo discharge of cohesive powders[J]. Powder Technology, 2009, 191(3): 272－279.

[28] 韓高孝，宮全美，周順華. 摩擦型巖土材料土拱效應(yīng)微觀機(jī)制顆粒流模擬分析[J]. 巖土力學(xué)，2013，34(6)：1791－1798.

Han Gaoxiao, Gong Quanmei, Zhou shunhua. Analysis of microcosmic mechanism of soil arching in frictional geotechnical material by particle flow simulating[J]. Rock and Soil Mechanics, 2013, 34(6): 1791－1798. (in Chinese with English abstract)

[29] Christoffersen J, Mehrabadi M M, Nemat-Nasser S. A micromechanical description of granular material behavior[J]. J Appl Mech, 1981, 48: 339－344.

[30] Fortin J, Millet O, de Saxcé G. Construction of an averaged stress tensor for granular medium[J]. European Journal of Mechanics A/Solids, 2003, 22: 567－582.

[31] 賴漢江，鄭俊杰，章榮軍，等．樁承式路堤土拱形成及荷載傳遞機(jī)制離散元分析[J]．巖土力學(xué)，2015，36(增刊1)：646－650.

Lai Hanjiang, Zheng Junjie, Zhang Rongjun, et al. Discrete element analysis of development and load-transfer mechanism of soil arching within piled embankment[J]. Rock and Soil Mechanics, 2015, 36(Supp.1): 646－650. (in Chinese with English abstract)

[32] 糧食鋼板筒倉設(shè)計規(guī)范GB50322-2011[S]. 北京：中國計劃出版社，2001.

[33] Hala Hammadeh. Flow and Pressure during Centric and Eccentric Discharge in Cylindrical Silos[J]. Science & Technology, 2012: 44－48.

[34] Koby?ka R, Molenda M. DEM modelling of silo load asymmetry due to eccentric filling and discharge[J]. Powder Technology, 2013, 233: 65－71.

Development of arching and silo wall pressure distribution in storage and discharging state based on discrete element analysis

Liu Kejin, Xiao Zhaoran, Wang Shihao

(450001)

It has been widely recognized that the most dangerous condition is the phenomenon that the emptying pressure on the wall of silos vibrates and increases obviously during discharging. However, the cause of this phenomenon is currently largely unknown. Most of the available papers and reports are based on continuum mechanics, that is, the stored material is regarded as a continuous entity on macroscopic level, and its particulate property is ignored. In fact, it is the microcosmic mechanical behavior of individual particles and the interaction between the particles and the silo walls that determine the emptying pressure distribution on the silo walls. Therefore, in this paper, the characteristics of emptying pressure were studied from the point of view of granular materials and microcosmic particle mechanics. A new method combining particle mechanics with discrete element method (DEM) was introduced to explore how the mechanical behavior of particles effects the distribution of emptying pressure. Firstly, the behavior of granular material in stored state was studied by DEM method. The silo with flat bottom is 0.5 m in diameter, 1.0 m in height and 0.1 m in outlet diameter, which is filled with 20 400 spherical particles. The distribution of the static wall pressure in stored state was verified by the test results and the Jassen Formula which is widely used in engineering. Secondly, the pressure distribution on the silo bottom wall was studied by simulating the discharge process, and the statistical analysis of multiply simulation results was performed, the arching effect near the outlet were proved according to the periodic pressure profile. Thirdly, in order to study the arching effect, three time points, i.e. start of arching, completion of arching and arch collapse, were selected on the periodic pressure profile. For each time point, the behaviors, such as force chain network, vertical stress distribution, and lateral stress distribution, the direction of principal stress, velocity field of granular material and so on, were systematically studied. On this basis, the features and evolution mechanism of the arching effect were investigated from the viewpoint of particle mechanics. Finally, the static pressure under the static stored state and the emptying pressure in discharging process on the silo wall were analyzed and verified by a model test, Standards (GB50322-2011) and published results. The results indicated that the arch which is 4.0 times wide of outlet diameter and 2.5 times high of outlet diameter was produced during the discharge process. Due to the formation of arching, the vertical stress above arching was transformed into the horizontal stress in a certain range above the arch foot, thus the horizontal stress was transferred to both sides of the silo wall and to the pressure of the silo wall was increased. When the arching collapsed, the vertical and horizontal above arching decreased, then the horizontal stress transferred to both sides of the silo wall also decreased, which resulted in to the decrease of wall pressure. The coefficient of overpressure reached a peak value of 2.70 at depth ratio of 0.35, while the peak wall pressure was 3.57 kPa at depth ratio of 0.85. It was found that the arching effect near the outlet behaved dynamically, following the rule of “start of arching-completion of arching-collapse of arch”, thus the dynamic arching effect is put forward. In addition, the important connection between the dynamic arching effect and the resulting emptying pressure distribution was identified. That is, the increased emptying pressure on silo wall is caused by the formation of arch, while the vibration of emptying pressure on wall of silo is induced by the dynamic arching effect. The research provides a new method of exploring the relationship between the arching effect and the emptying pressure distribution on wall of silos from the micro and macro aspects. The findings obtained in this paper can provide references for revealing the load-transfer mechanism from particles to silo wall.

stresses; numerical methods; silos; wall pressure; arch; soybean; granular materials; micromechanics

10.11975/j.issn.1002-6819.2018.20.035

S126

A

1002-6819(2018)-20-0277-09

2017-12-08

2018-07-19

國家自然科學(xué)基金項目（5170081095）；河南工業(yè)大學(xué)高層次人才基金（31400375）

劉克瑾，講師，博士。主要從事糧食顆粒的力學(xué)性質(zhì)、筒倉結(jié)構(gòu)設(shè)計和地基處理研究。Email：Kejin_Liu@126.com

劉克瑾,肖昭然,王世豪. 基于離散元模擬筒倉貯料卸料成拱過程及筒倉壁壓力分布[J]. 農(nóng)業(yè)工程學(xué)報，2018，34(20)：277－285. doi：10.11975/j.issn.1002-6819.2018.20.035 http://www.tcsae.org

Liu Kejin, Xiao Zhaoran, Wang Shihao. Development of arching and silo wall pressure distribution in storage and discharging state based on discrete element analysis[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(20): 277－285. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2018.20.035 http://www.tcsae.org