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

劉克瑾，肖昭然，王世豪

?

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

劉克瑾，肖昭然，王世豪

（河南工業大學土木建筑學院，鄭州 450001）

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

應力；數值方法；筒倉；倉壁壓力；拱；大豆；顆粒物質；細觀研究

0 引 言

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

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

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

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

1 模型建立

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

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

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

表1 筒倉離散元法模型材料參數

2 貯料工況分析

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

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

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

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

圖3 貯料靜工況顆粒物質的力鏈網絡、孔隙率、豎向應力和橫向應力分布

對比圖3a和3b發現，在顆粒塊體內部，力鏈強度較弱且分布密集。弱力鏈在整個筒倉內都有分布，而強力鏈主要分布在孔隙率較小的區域。對比圖3a、3b和3c發現，在豎向應力較大的區域內，強度較大的力鏈方向與重力方向夾角較小；在橫向應力較大的區域內，強度較大的力鏈方向與重力方向夾角較大。說明強度較大的力鏈的方向影響區域內豎向應力及橫向應力的分布。

3 卸料過程

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

3.1 倉底壓力分析

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

圖4 拱腳位置

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

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

1）卸料過程中倉底及倉壁的壓力一直處于震蕩變化狀態，這一現象與文獻[13,25-26]的研究相符。靠近卸料口的倉底壓力較小，5號測墻倉底壓力與倉底峰值壓力之比不超過0.41，4號測墻不超過0.39。遠離卸料口倉底壓力較大，如1、2號和7、8號測墻，倉底壓力與倉底峰值壓力之比在0.5與1.0之間周期性變化。

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

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

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

3.2 細觀顆粒層次的拱效應分析

3.2.1 結拱起始

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

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

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

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

分析可知：在強度較大的力鏈分布的區域內，豎向應力或橫向應力較大，這是因為應力的分布是由區域內力鏈強度決定的。對比圖6b、6c和6d發現，橫向應力的分布與區域內顆粒速度無明顯聯系；而在豎向應力較大的區域內，顆粒速度較小，反之亦然，說明區域內顆粒速度主要由豎向應力的分布決定。

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

3.2.2 結拱完成

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

觀察圖7可以發現：

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

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

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

圖7 結拱完成時顆粒物質的力鏈網絡、豎向應力、橫向應力和速度分布

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

3.2.3 拱塌落

圖8為拱塌落時刻點C的貯料力鏈網絡分布圖，豎向應力分布圖，橫向應力分布圖及顆粒速度圖。

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

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

圖8 拱塌落時顆粒物質的力鏈網絡、豎向應力、橫向應力和速度分布

3.3 主應力方向分析

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

圖9 不同時刻主應力方向分布圖

圖9a顯示，貯料工況下主應力方向以豎向為主，且除筒倉頂部主應力方向較雜亂外，其他區域內主應力方向分布有較好的規律性：距倉壁越近，主應力方向與重力方向夾角越大。圖9b顯示，結拱起始時刻，筒倉下部側壁附近主應力方向以豎向為主，但其他區域內的主應力方向較貯料工況發生了較大偏轉。圖9c顯示，結拱完成時刻，在懸鏈線處形成了弧形拱狀結構，拱下方的主應力方向分布較為雜亂，在拱頂的正上方，由于拱的形成阻礙了上方顆粒的流動，對上方顆粒起到了支撐作用，拱頂上方的主應力方向以豎向為主。圖9d拱塌落時刻，除筒倉下部側壁附近主應力方向以豎向為主，有著較好的規律性外，其他區域內主應力方向分布較為雜亂。

對比圖9b、9c和9d發現，雖然圖9b和9d中懸鏈線內部主應力方向有拱形分布，但其附近的主應力方向分布較為雜亂，說明懸鏈線附近沒有拱的形成，而圖9c中懸鏈線附近的主應力方向分布較圖9b和9d有更好的規律性，說明懸鏈線附近有拱的形成。在拱的形成過程中，由于拱對上方顆粒的支撐作用，拱頂上方的主應力方向與重力方向的夾角減小；在拱的塌落過程中，由于失去了拱的支撐作用，拱頂上方的主應力方向發生偏轉，與重力方向的夾角增大。

3.4 宏觀倉壁卸料壓力分析

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

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

圖10 倉壁側壓力分布

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

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

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

3.5 細觀動態成拱機制對宏觀倉壁卸料壓力分布的影響

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

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

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

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

4 結 論

本文采用離散元方法和室內模型試驗模擬了筒倉貯料工況及卸料工況，模型試驗結果驗證了離散元法結果。通過系統的研究結拱起始、結拱完成和拱塌落等時間節點的貯料力鏈網絡、豎向應力和橫向應力、顆粒速度及倉壁壓力分布，從細觀顆粒層次探討了卸料過程中的動態成拱機制，以及由此引發的宏觀倉壁側壓力分布。主要得出以下結論：

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

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

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

4）拱的形成導致卸料過程中倉壁壓力增大，拱塌落使得倉壁壓力減小。動態成拱機制即“結拱（結拱起始-結拱完成）-拱塌落”導致倉壁的卸料壓力呈現“較小-局部極大-較小”的震蕩分布。

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

6）本文的動態壓力修正系數最大值為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] 劉定華，郭明昌，楊建斌. 筒倉倉壁動態壓力的測試和分析[J]. 西安建筑科技大學學報：自然科學版，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] 李國柱，劉定華. 筒倉動態壓力的計算和測試[J]. 寧波高等專科學校學報，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] 趙松. 筒倉貯料壓力分析及其應用[D]. 武漢：武漢理工大學，2013.

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

[7] 車穎文. 筒倉的受力情況和穩定性研究[D]. 武漢：武漢理工大學，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] 孫其誠，王光謙. 靜態堆積顆粒中的力鏈分布[J]. 物理學報，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] 鄒同彬，陸渝生，連志穎. 動光彈實驗中的主應力分離方法[J]. 解放軍理工大學學報：自然科學版，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] 譚援強，肖湘武，鄭軍輝，等. 錐形改流體下部孔徑對筒倉卸料流態的影響[J]. 農業工程學報，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]. 農業工程學報，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] 肖昭然，王軍，何迎春. 筒倉側壓力的離散元數值模擬[J]. 河南工業大學學報：自然科學版，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] 張昭，劉克瑾，肖昭然，等. 筒倉倉壁摩擦對倉壁側壓力影響的研究[J]. 河南工業大學學報：自然科學版，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] 馮家暢，程緒鐸，嚴曉婕. 大豆與倉壁材料摩擦系數的研究[J]. 大豆科學，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] 程緒鐸，陸琳琳，石翠霞，等. 大豆內摩擦角的測定與實驗研究[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]. 中國糧油學報，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] 孫其誠，辛海麗，劉建國，等. 顆粒體系中的骨架及力鏈網絡[J]. 巖土力學，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] 房營光，侯明勛，谷任國，等. 樁承式路堤中土拱效應產生過程可視化分析[J]. 巖土工程學報，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]. 土木工程與管理學報，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] 張大英，許啟鏗，王樹明，等. 筒倉動態卸料過程側壓力模擬與驗證[J]. 農業工程學報，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] 韓高孝，宮全美，周順華. 摩擦型巖土材料土拱效應微觀機制顆粒流模擬分析[J]. 巖土力學，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]．巖土力學，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] 糧食鋼板筒倉設計規范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

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

劉克瑾，講師，博士。主要從事糧食顆粒的力學性質、筒倉結構設計和地基處理研究。Email：Kejin_Liu@126.com

劉克瑾,肖昭然,王世豪. 基于離散元模擬筒倉貯料卸料成拱過程及筒倉壁壓力分布[J]. 農業工程學報，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