深淺倉卸料壓力離散元數值模擬研究

原 方，劉海林，程遠浩，杜 乾

1.河南工業大學 土木建筑學院，河南 鄭州 450001 2.三明學院 建筑工程學院，工程材料與結構加固福建省高等學校重點實驗室，福建 三明 365004

糧食儲備庫中筒倉儲糧占絕大多數，不僅如此，在煤炭、冶金、電力、建材等行業，越來越多的筒倉也被興建，因此，筒倉安全的重要性不言而喻。

筒倉卸料壓力增大機理復雜多變，一直困擾著國內外學者。許多國家的筒倉規范中計算動態側壓力的方法幾乎全部采用靜壓力乘以一個側壓力的修正系數[1-2]，具體方法為：關于深倉動壓力，筒倉設計者大多采用Janssen理論，但設計使用的修正系數不同；關于淺倉動壓力，多采用Rankine理論和修正后的Coulomb理論。側壓力修正系數過大會造成工程材料的浪費，過小又會威脅筒倉結構安全，因此用合適的側壓力修正系數來準確預測動壓力亟須被解決。針對這一難題，劉定華等[3]提出了由整體流動引起筒倉側壁動壓力的計算方法；原方等[4]提出了考慮倉徑比影響的一種新的淺圓倉散料側壓力的計算方法，已被筒倉規范采用，主要來解決大型淺圓倉側壓力的計算;孫珊珊等[5]提出了一個基于統一強度理論，考慮3個主應力的共同影響，深倉和淺倉都適用的統一側壓力系數。此外關于超壓系數(卸料時的動壓力與儲料時的靜壓力之比)，研究發現，卸料過程中動壓力達到靜壓力的2～3倍[6-7]，這是流態的變化或動力拱的形成引起的[8-12]。但截止到目前，尚未出現被廣泛接受的合理的理論來準確計算筒倉側壁上的動壓力，仍未搞清動壓力的增大機理，無法定性定量的解決這一科學問題。

作者在試驗驗證的基礎上，運用離散元軟件PFC2D建立6個高徑比(H/D)不同的模型筒倉進行數值模擬，其中3個深倉(H/D>1.5)高徑比分別為2.9(模型倉a)、2.2(模型倉b)和1.76(模型倉c),3個淺倉(H/D<1.5)高徑比分別為1.47(模型倉d)、1.28(模型倉e)和1.1(模型倉f)，對靜壓力、動壓力以及超壓系數進行對比分析，另外從流態和動力拱兩個不同的角度來分析筒倉側壁動壓力增大的情況，為能準確預測筒倉側壁動壓力的計算公式提供推導依據，為構建筒倉結構安全設計的機理研究提供參考。

1 模型倉試驗和數值模擬

1.1 試驗模型倉

試驗采用有機玻璃筒倉，便于觀察卸料過程中小麥的流態，底部漏斗高0.7 m，筒倉倉壁高5.5 m，裝料高4.4 m。配套設備有中轉倉、斗提機和傳送帶，可以實現玻璃筒倉入料和卸料的機械化，如圖1所示。試驗所用小麥重力密度為8 kN/m3，內摩擦角為25°，外摩擦系數為0.4，小麥物理參數選取參考文獻[8]。

1. 筒倉筒體; 2. 帶孔漏斗; 3. 鋼制筒臺; 4. 平臺鋼架圖1 模型筒倉及結構圖Fig.1 Model silo and structure

試驗采用中心卸料方式，傳感器布置如圖2所示。A、B、C每列分別布置10個土壓力傳感器，D列布置4個薄膜傳感器。

圖2 倉壁傳感器立面圖Fig.2 Sensor elevation layout of silo wall

1.2 模擬模型倉

模擬模型倉以試驗筒倉為原型，裝料高度4.4 m，模型調試采用控制變量法，通過對模型試驗數據進行反復標定以確定模型的各項物理參數[13]，模擬模型倉的各項主要材料參數為：墻體法向剛度1.0×107Pa，墻體切向剛度0.5×107Pa，顆粒法向剛度1.0×106Pa，顆粒切向剛度0.5×106Pa，顆粒重力密度800 kg/m3，顆粒外摩擦系數0.4，顆粒內摩擦系數0.5，模型倉a測點深度及滿倉狀態如圖3所示。

圖3 模型倉a測點深度及滿倉狀態Fig.3 Depth and full filling state of tested point ain a model silo

1.3 試驗結果與模型倉a模擬結果對比分析

圖4中試驗靜壓力和動壓力是A、B和C 列中傳感器所得平均值，由圖4可知，模擬的靜壓力與試驗的靜壓力總體趨勢一致。試驗最大動壓力發生在筒倉下部大約1/3處，為12.23 kPa，模擬最大動壓力發生在筒倉下部大約1/4處，為18.26 kPa。試驗超壓系數是A、B和C 列中傳感器平均動壓力與平均靜壓力的比值,結果如圖5所示，試驗最大超壓系數發生在筒倉中下部位置，為2.05，模擬最大超壓系數發生在筒倉下部位置，為2.48。由圖4和圖5可知，試驗結果與模擬結果不盡相同，一方面是由于計算機技術的限制，小麥粒徑設置過小則無法進行計算，數值模擬時所采用的粒徑為8 mm，與試驗倉所使用的小麥粒徑不同；另一方面是由于模擬模型倉側壁與漏斗相交處附近，邊界條件設置復雜，文中為簡化模型，此處設置為線性接觸，但誤差在可接受的范圍內。

圖4 試驗與模擬的靜壓力、動壓力對比Fig.4 Comparison of static pressure and dynamic pressure between experiment and simulation

圖5 試驗超壓系數與模擬超壓系數對比Fig.5 Comparison of experimental overpressure coefficient and simulated overpressure coefficient

2 高徑比不同的深淺倉數值模擬

2.1 高徑比不同的深淺倉模擬模型

所有的模擬模型倉裝料高度、測點布置深度以及建模時所用的材料屬性參數相同，不同的只是筒倉直徑以及卸料口寬度。卸料口寬度為筒倉直徑的1/5，底部漏斗傾角50°，模型筒倉裝糧高度4.4 m，模型倉a倉徑1.5 m，模型倉b倉徑2 m，模型倉c倉徑2.5 m，模型倉d倉徑3 m，模型倉e倉徑3.5 m，模型倉f倉徑4 m。

2.2 高徑比不同的深淺倉模擬結果分析

從PFC2D軟件中提取模型倉b卸料壓力曲線如圖6所示，橫坐標表示筒倉卸料時步(時步為PFC2D軟件中進行迭代得到最終狀態的計算單位)，前期曲線的平線段表示各個不同深度測點的靜壓力，由圖6可知，最大卸料壓力發生在卸料初期，大約在2.5×105時步處，發生超壓且現象明顯，且在卸料過程中，同一高度的卸料壓力并不是一直增大或減小，而是出現上下震蕩。其他模型倉與之類同。

圖6 PFC中模型倉b卸料壓力曲線Fig.6 Discharge pressure curve of model b in PFC

圖7 高徑比不同的5個模型倉靜態壓力值對比Fig.7 Comparison of static pressure values of five model silos with different height diameter ratios

圖8 高徑比不同的5個模型倉動態壓力值對比Fig.8 Comparison of dynamic pressure values of five model silos with different height diameter ratios

圖9 高徑比不同的5個模型倉超壓系數對比Fig.9 Comparison of overpressure coefficient of five model silos with different height-diameter ratios

高徑比不同的5個模型倉靜壓力對比、動壓力對比、超壓系數對比如圖7—圖9所示，可以發現：(1)動壓力明顯大于靜壓力，筒倉下部位置最為明顯。(2)深倉模型與淺倉模型動壓力相差不大，模型倉b最大動壓力16.42 kPa，模型倉c最大動壓力18.43 kPa，模型倉d最大動壓力19.21 kPa，模型倉e最大動壓力19.01 kPa，模型倉f最大動壓力18.03 kPa，由前文知模型倉a的最大動壓力18.36 kPa，且最大動壓力均發生在筒倉下部位置。可以看出，相同高度、較大直徑的筒倉最大動壓力隨直徑增大呈整體上升趨勢，但非線性，這與付建寶[14]發現的一致，同時也說明了設計淺倉時，只采用Rankine理論和修正后的Coulomb理論作為設計值，不考慮動態壓力的修正是不科學的。

3 動態壓力增大原因分析

3.1 從貯料流態的角度分析

顆粒在離散元軟件PFC2D中分層生成，結果如圖10所示。由于模擬數據較多，故下文中卸料過程、速度場和力鏈場，深倉以模型倉a為例，淺倉以模型倉f為例。

注：(a)、(b)、(c)、(d)表示時步分別為2.3×105、3.3×105、5.3×105、6.3×105，圖11—圖16同。圖10 筒倉內散體顆粒生成過程Fig.10 Formation process of bulk particles in silos

筒倉內貯料在重力作用下流動，可分為整體流動和管狀流動，整體流動為筒倉內貯料均勻向下運動，同一高度的貯料保持在同一水平上；由于筒倉下部截面形態的變化阻礙了筒壁處貯料的運動，導致筒壁處的貯料流動緩慢或存在短暫的死料區，使得筒倉下部處于同一水平面上的貯料出現逐漸增大的速度梯度，從而形成管狀流動。如圖11所示，模型倉a卸料前中期，筒倉內貯料呈整體下落態勢，即每層顆粒基本上保持矩形形狀不破壞，在重力作用下整體下落；筒倉內貯料大概剩余初始1/3時，開始出現管狀流動，中間顆粒流動速度明顯大于兩側，上表面開始凹陷，卸料逐漸從單一的整體流動變為整體流動與管狀流動混合的卸料方式，且與試驗觀察到的流態高度一致。如圖12所示，模型倉f卸料前期，筒倉內顆粒整體流動，筒倉內貯料大概剩余初始2/3時，開始出現管狀流動，與深倉模型a一致，開始出現混合流，并且可以發現，當筒徑足夠大時，整體流動狀態不再出現。整體流動是筒倉設計首選流態，具有料流穩定、不結塊以及不分離的優點，通過改流裝置使管狀流動轉變為整體流動，且發現筒倉側壁壓力峰值減小[15]，故前文所述淺倉相對于深倉卸料壓力峰值偶有增大，兩者流態的不同是卸料壓力不同的重要原因之一。但Mathews等[16]也指出很難在整個卸料過程實現“先進先出”的整體流，因倉壁與貯料的剪切力極易造成貯料流動的速度梯度，從而在卸料末期形成漏斗流。

圖11 深倉模型a卸料過程Fig.11 Discharging process of deep silo model a

圖12 淺倉模型f卸料過程Fig.12 Discharging process of shallow silo model f

如圖13和圖14所示，無論模型倉a還是模型倉f，隨著卸料時間的延長，筒倉中間的顆粒流動速度要比兩側顆粒流動速度快，導致模型倉a卸料后期，模型倉f中后期上層貯料最先流出，這與上文所討論的流態相一致。

圖13 深倉模型a卸料速度場Fig.13 Discharging velocity field of deep silo model a

圖14 淺倉模型f卸料速度場Fig.14 Discharging velocity field of shallow silo model f

3.2 從貯料結拱的角度分析

筒倉內貯料在倉壁截面變化處形成擠壓流動，使貯料在該處不斷形成結拱-塌拱的狀態，即不斷形成動力拱，從圖15和圖16可以看出，無論是深倉還是淺倉，筒倉側壁附近的力鏈隨著筒倉深度的加深而變得密集，即接觸力越來越大，并且在卸料過程中，筒倉中心的貯料流動速度快但接觸力小，貯料迅速減少形成的動力拱則會抑制兩側的貯料，兩側的貯料流動速度慢但接觸力大，這種速度的不同造成了力學梯度，使筒倉內的貯料在流動時受到剪切混合力的作用，分布形態成非水平狀態即力鏈在水平方向上分布的不均勻, 是造成壓力凹陷[17]的主要原因。與王學文[9]發現的情況一致，并且結拱所需時間與筒倉尺寸有關，對于圓筒倉而言，直徑越大貯料顆粒在卸料過程中結拱所需時間也就越長，貯料顆粒本身還未來得及調整結拱的最佳位置卻已經下移，那么結拱的可能性就越小,這是卸料時淺倉相對于深倉卸料壓力不同的重要原因之一。

圖15 深倉模型a卸料力鏈場Fig.15 Discharging force chain field of deep silo model a

圖16 淺倉模型f卸料力鏈場Fig.16 Discharging of focus chain fieldshallow silo model f

4 結論

從筒倉高徑比入手，采用室內模型試驗法和離散單元法研究貯料在靜態儲糧狀態和卸料過程中的力學行為。試驗最大動壓力發生在筒倉下部約1/3處，為12.23 kPa，模擬最大動壓力發生在筒倉下部約1/4處，為18.36 kPa；試驗最大超壓系數發生在筒倉中下部，為2.05，模擬最大超壓系數發生在筒倉下部，為2.48。經過對高徑比不同的6個深淺倉模型的結果分析，相同高度、較大直徑的筒倉最大動壓力隨直徑增大呈整體上升趨勢，但非線性。在卸料過程中，筒倉最大卸料壓力發生在卸料初期，同一高度的卸料壓力并不是一直增大或減小，而是在卸料壓力的增大或減少的趨勢下出現上下震蕩。深倉模型a內貯料大概剩余初始1/3時，淺倉模型f內貯料大概剩余2/3時，開始出現管狀流動，中間顆粒流動速度明顯大于兩側，上表面開始凹陷，卸料逐漸從單一的整體流動演變為整體流動與管狀流動混合的卸料方式，兩者流態的不同是卸料壓力不同的重要原因之一。無論是深倉還是淺倉，筒倉側壁附近的力鏈隨著筒倉深度的加深而變得密集，即接觸力越來越大，并且在卸料過程中，筒倉中心的貯料流動速度快但接觸力小，貯料迅速減少形成的動力拱則會抑制兩側的貯料，兩側的貯料流動速度慢但接觸力大，并且結拱所需時間與筒倉尺寸有關，對于圓筒倉而言，直徑越大貯料顆粒在卸料過程中結拱所需時間也就越長，貯料顆粒本身還未來得及調整結拱的最佳位置卻已經下移，那么結拱的可能性就越小,這是卸料時淺倉相對于深倉卸料壓力不同的重要原因之一。