成拱工況下倉壁超壓彈性波動態(tài)模擬分析

馮 永，原子然，張盼盼，謝飛亞，李曉一，李 萌

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

近年來，卸糧成拱對倉壁造成破壞的現(xiàn)象時有發(fā)生，但現(xiàn)有室內實驗及理論計算模型無法客觀反映筒倉卸糧成拱下的倉壁超壓機制形成過程。為了適應糧食行業(yè)的迅速發(fā)展，深入研究倉壁側壓力勢在必行。

國內外學者提出了理論計算公式量化筒倉倉壁壓力，如Jassen公式，但并不適用于卸糧工況，更不適用于成拱工況[1-10]。部分學者研究了成拱時側壁產生的超壓，發(fā)現(xiàn)物料含水率、物料的內外摩擦系數(shù)、料倉尺寸等均會影響卸料成拱[11-14]。卸料過程中的成拱現(xiàn)象會導致側壁超壓[15]，拱內顆粒承受倉中粒料的大部分應力。實際工程中，倉壁并不是絕對剛性的，倉壁彈性變形傳播在卸料過程中影響倉壁動壓力的形成[16-18]，但是超壓形成的細觀機制尚不明確。

在卸糧成拱這一工況中，研究倉壁壓力彈性波分布規(guī)律對于揭示倉壁超壓作用機理具有重要借鑒意義。本文中，通過工程實例建立縮尺模型，運用室內實驗和PFC軟件對筒倉卸料過程進行動態(tài)模擬，揭示筒倉卸糧成拱工況下超壓彈性波的形成對倉壁超壓影響，分析二次超壓的形成和傳播機制，為預防類似事故提供理論依據。

1 形成機理

1.1 成拱實現(xiàn)過程

糧食是一種散粒體，散粒體兼具液體和固體的一部分特性，但是其物理力學特性與固體液體又都不盡相同。從細觀的角度來看，散體的本質是很多離散顆粒的集合，在一定范圍內呈現(xiàn)出與液體相類似的流動性，從整體上來看依然可以認為貯料是各向同性的散粒集合體。

貯料在靜態(tài)工況下是不可能成拱的，成拱的過程其實是一個動態(tài)平衡的過程。在卸料過程中，倉內顆粒在自重作用下開始下瀉，同時，顆粒之間以及顆粒和倉壁之間的摩擦擠壓碰撞等耦合作用使得顆粒下瀉呈現(xiàn)復雜的運動狀態(tài)，顆粒的位置和速度以及所受到的合力和合力矩每時每刻都在變化。各顆粒之間的相互移動和摩擦擠壓會使顆粒通過位置和應力狀態(tài)的不斷調整找到一個可以互相傳遞壓力的最佳位置，形成連接倉壁兩側的完整力鏈，且力鏈之間不產生顆粒的空位，以實現(xiàn)力向拱腳的傳遞和對上部物料流動的阻礙，并要求拱腳有足夠的強度用于支撐。由于力鏈的形成影響了顆粒的正常出流，因此顆粒流動速度的降低產生與顆粒流速方向一致的慣性力，筒倉內拱形力鏈以上的顆粒自重以及新產生的慣性力全部作用于顆粒拱，并通過拱腳傳遞到倉壁兩側，使得拱腳處倉壁產生超壓。

1.2 超壓彈性波形成機理

在實際工程和數(shù)值模擬中，倉壁都不是理想的絕對剛性體，倉壁具有一定的彈性，拱腳處的超壓會使緊貼拱腳的倉壁產生微小的橫向膨脹變形，且超壓會造成儲料密度增大，膨脹段容積的增大和儲料擠密造成的糧食體積縮小會容納更多的貯料，流入的儲料也具有一定速度，在流入膨脹區(qū)段后必然伴隨速度降低為零的過程和隨之產生的慣性力。倉壁膨脹示意圖如圖1所示，隨著貯料的不斷流入擴張體積中，成拱超壓自下而上以波的形式在倉壁上傳播。

圖1 倉壁膨脹示意圖Fig.1 Diagram of silo wall expansion

2 工程實驗相似性認證

根據工程實例[11]建立實驗模型，實驗與模擬均采用20∶1縮尺模型，參與卸糧成拱過程的主要物理參數(shù)有重力加速度g、糧食密度ρ、筒倉模型幾何尺寸L、墻體彈性模量E、糧食黏聚力C、糧食內部的摩擦力φ1、糧食外部的摩擦力φ2、倉壁受到的應力為σ。其一般函數(shù)式為

通過相似性準則可以得到函數(shù)表達式

幾何相似性系數(shù)CL為1/20；因為密度和重力加速度難以在實驗中改變，所以采用密度相似系數(shù)Cρ=1、重力加速度相似系數(shù)Cg=1。根據Jassen公式估算應力相似性系數(shù)，并帶入式（2）中相似性準則驗證，可以認證本文中模型試驗的材料相似性、幾何相似性與應力相似性之間相互滿足，從而證明所采用的縮尺模型滿足相似性認證。

3 實驗研究

3.1 實驗模型及壓力波監(jiān)測

該工程實例中破壞倉的內徑為Φ6 000 mm，倉壁厚160 mm，倉壁高度為20 000 mm，卸料口大小為Φ1 000 mm。成拱實驗模型結構示意圖如圖2所示。

圖2 成拱實驗模型結構示意圖Fig.2 Schematic diagram of arch-forming experimental model

成拱實驗所用20∶1的縮尺模型尺寸如表1所示。

表1 成拱實驗模型尺寸Tab.1 Size of the arch model

筒倉用厚度為8 mm的亞力克材料制作，以期清楚地觀測到顆粒的流動狀態(tài)。將實驗模型倉內裝滿粒徑為5～6 mm的小麥，待裝料完成后撤掉料斗底板進行卸糧實驗，通過添加粉塵、改變含水率、倉壁貼砂紙等方法，進行多次實驗，最后實現(xiàn)筒倉卸糧成拱，成拱實驗系統(tǒng)如圖3所示。

圖3 成拱實驗系統(tǒng)Fig.3 Experimental system of arched grain

在筒倉左側倉壁布置10個監(jiān)測點，每個監(jiān)測點上安放一個應變片來監(jiān)測倉壁側壓力，可以得到倉壁側壓力在成拱過程中的變化圖像，倉壁監(jiān)測點分布如圖4所示。

圖4 倉壁監(jiān)測點分布圖Fig.4 Distribution of monitoring points on silo wall

3.2 實驗結果分析

實驗發(fā)現(xiàn)，成拱現(xiàn)象大致發(fā)生在筒倉的1/3高度處，1/3高度以下的顆粒自由卸料，而筒倉上部在卸料過程很快迎來了零壓力區(qū)，取0.1 m2墻段只繪制編號為 4、5、6、7、8 這 5 個監(jiān)測點側壓力圖線，各監(jiān)測點處側壓力如圖5所示，可以看到，成拱過程中側壁成拱壓力的形成過程。伴隨著卸糧成拱，倉側壓力波動，并在波峰出現(xiàn)穩(wěn)定側壓，即為成拱壓力，但受限于客觀實驗條件不能精確觀測倉側成拱壓力形成的細觀機理，因而在后文中通過PFC仿真模擬研究。

圖5 倉壁監(jiān)測點處側壓力Fig.5 Pressure at monitoring points

4 數(shù)值分析模型

4.1 線性接觸模型與本構模型

4.1.1 線性接觸模型

球與球之間以及球與墻面之間的接觸都采用線性接觸模型。線性接觸規(guī)定了兩物體之間無窮小的交界面是不能抵抗彎矩的，接觸彎矩常等于0。線性接觸中由線性元件和阻尼元件來提供接觸力，線性元件提供無張力的線彈性力和摩擦力，而阻尼元件提供2個單元體之間的粘結力，將線性元件基本等同于一個彈簧和滑移板，將阻尼元件等同于一個阻尼器，線性接觸模型如圖6所示。

圖6 線性接觸模型Fig.6 Linear contact model

4.1.2 本構模型

1）顆粒與顆粒

根據力-位移定律，得到2個顆粒接觸時的單位矢量模型和顆粒接觸變形的情況。2個糧食顆粒i、j相互接觸的示意如圖7所示，設2個顆粒中心坐標分別為（xi，yi）、（xj，yj），因為小麥顆粒的粒徑分布雖然不同，但是差異不大，可以認為2個顆粒半徑均為 rc，接觸點處的法向向量n→=（nx，ny），切向向量為t→=（tx，ty）。

圖7 糧食顆粒接觸示意圖Fig.7 Schematic diagram of contact between grain particles

由此可以得接觸點處的法向和切向向量的表達式為

而顆粒i、j接觸點的速度分別表示為

將式（6）、（7）帶入式（5），可得到

相對速度在法向和切向的分量為

對于顆粒i、j，法向切向的速度分量都造成了相應方向的位移，根據上述方式所計算出的速度分量，可以得到一個時間段變化的位移分量的表達式為

根據力-位移定律可以得到力位移關系的表達式為

由此易得顆粒i的運動方程為

式中I為顆粒的慣性矩，在線性接觸模型中一般認為顆粒接觸點的彎矩常等于零。

在運算的每個微小時段，細觀參數(shù)都在發(fā)生改變。以顆粒i為例，力在一個微小時段的更新為

速度在一個微小時段的更新為

顆粒位置在一個微小時段的更新為

2）顆粒與墻體

顆粒與墻體的接觸示意圖如圖8所示。設代號為i的顆粒中心點坐標為（xi，yi），代號j的墻體圍繞點（xj，yj）發(fā)生旋轉。墻壁與顆粒接觸點坐標為（xc，yc）。接觸點處的法向向量=（nx，ny），切向向量為=（tx，ty）。

圖8 顆粒與墻體接觸示意圖Fig.8 Schematic diagram of contact between particle and wall

而顆粒i與墻體j在接觸點的速度分別為

式中，假定墻體不發(fā)生平動，所以vj=0為球體的半徑向量，為墻體轉動的半徑向量，wj為墻體轉動角速度。

將式（23）、（24）帶入式（22），可以得到接觸點處的相對移動速度

相對速度在法向和切向的分量為

接觸力法向切向的增量及位置速度和力的更新方法均依前文所述，不再贅述。

4.2 數(shù)值模型的建立

參照模型實驗建造同樣是20∶1的縮尺模型，根據文獻[12]所提供的糧食力學特性數(shù)據，并參照上文所述的相似條件，選取利于成拱的顆粒和墻壁的力學參數(shù)，如表2所示。Wall和Ball分別表示墻體單元和顆粒單元，kn、ks則分別為它們的法向和切向剛度。

表2 模擬中顆粒和墻壁力學參數(shù)Tab.2 Mechanical parameters of particle and wall in simulation

墻體和顆粒的線性接觸模型所選取的參數(shù)如表3所示。

表3 模擬中接觸模型的主要參數(shù)Tab.3 Main parameters of the contact model in simuloction

在PFC模型中，倉壁由分段建立的20個Wall單元組成，分別代表20個虛擬監(jiān)測點監(jiān)測成拱工況下倉壁動態(tài)側壓力。分層裝料完成并達到的穩(wěn)定裝載模型如圖9所示。

圖9 穩(wěn)定裝載模型Fig.9 Steady loading model

4.3 動態(tài)監(jiān)測點的布置

顆粒在筒倉卸糧成拱的過程中呈現(xiàn)復雜的，屬于非對稱、非均勻的流動狀態(tài)。參照上述試驗模型沿左右倉壁各生成10段100 mm的墻體，代表設置10個虛擬監(jiān)測點以監(jiān)測卸糧過程中倉壁壓力分布變化，仿真模擬中的監(jiān)測點布置如圖10所示。

圖10 仿真模擬中的監(jiān)測點布置Fig.10 Aarrangement of monitoring points in simulation

5 數(shù)值模擬分析與結果

通過采用有利于成拱的物理力學參數(shù)進行多次模擬試驗，最終實現(xiàn)成拱。對自重裝料并達到穩(wěn)定后進行的卸料過程進行動態(tài)模擬分析，卸糧成拱時的顆粒分布圖如圖11所示。

圖11 卸糧成拱時的顆粒分布圖Fig.11 Distribution of particles when grain is arched

5.1 力鏈形成過程

在卸料過程中，顆粒的速度位置和應力不斷發(fā)生變化，逐漸在筒倉內形成一個連接倉壁兩側的不規(guī)則拱形力鏈。如果顆粒之間可以很好地傳遞應力而且拱腳擁有足夠的強度，就會出現(xiàn)卸糧成拱現(xiàn)象，阻礙了拱形上方粒料的出流。成拱過程中的力鏈變化如圖12所示。

圖12 成拱過程中力鏈變化Fig.12 Force chain formation process

5.2 超壓彈性波形成的模擬分析

筒倉成拱的拱腳位置大約在距離筒倉和料斗交界約1/3處，因此在3號監(jiān)測點以上位置的側壓力受到筒倉成拱超壓的影響，而1、2、3號監(jiān)測點的側壓力變化均未受到成拱的影響；又因為成拱時已經卸出一部分顆粒，倉壁出現(xiàn)零壓力區(qū)，所以選取4、5、6、7、8共5個監(jiān)測點進行觀察和分析。

PFC軟件可以直接繪制出倉壁側壓力隨時步變化的圖像，但是無法精確看到每一時步所對應的側壓力數(shù)值，所以直接使用PFC輸出每一時步各點對應側壓力的Excel表格，并繪制出動態(tài)側壓力的折線圖。

左側倉壁各監(jiān)測點的動態(tài)側壓力如圖13所示。倉壁成拱左側壓力分布如圖14所示，成拱超壓系數(shù)如表4所示。

圖13 左側倉壁動態(tài)側壓力Fig.13 Ddynamic pressure at each monitoring points on left side of arch

圖14 倉壁成拱左側壓力分布圖Fig.14 Pressure distribution on left side of arch

表4 左側倉壁各監(jiān)測點側壓力超壓系數(shù)Tab.4 Overpressure coefficient at each monitoring point on left side of arch

仿真模擬中，右側倉壁各監(jiān)測點的動態(tài)側壓力如圖15所示。

圖15 右側倉壁動態(tài)側壓力Fig.15 Dynamic pressure at each monitoring points on right side of arch

倉壁成拱右側壓力分布如圖16所示，成拱超壓系數(shù)如表5所示?？梢郧逦乜吹?，倉壁側壓力的分布隨著深度的增加持續(xù)增大，而成拱時的超壓系數(shù)也隨深度而增加，且成拱時側壓力的超壓并不出現(xiàn)在同一時間點，這應該與在拱形上方由下而上速度降低的先后順序有關，也與成拱后倉壁由拱腳部分自下而上以彈性波的形式傳播超壓的行為有關。

圖16 倉壁成拱右側壓力分布圖Fig.16 Pressure distribution on right side of arch

表5 右側倉壁各監(jiān)測點側壓力超壓系數(shù)Tab.5 Overpressure coefficient at each monitoring point on right side of arch

對比圖5、13、15可以發(fā)現(xiàn)，側壓力隨卸料時間變化的趨勢與室內試驗結果相符合，實驗中成拱側壓力的大小與模擬結果相差不大，由此可以證明PFC模擬的精確性。

5.3 超壓彈性波傳播規(guī)律的模擬分析

模擬中，成拱超壓時左側倉壁各監(jiān)測點時步如圖17所示。

圖17 成拱超壓時左側倉壁各監(jiān)測點出現(xiàn)的時步Fig.17 Time step of arched overpressure occurrence at monitoring points on left side of arch

模擬得到的成拱超壓時右側倉壁各監(jiān)測點步如圖18所示。

圖18 成拱超壓時右側倉壁各監(jiān)測點的時步Fig.18 Time step of arched overpressure occurrence at monitoring points on right side of arch

鑒于筒倉結構的對稱性，多數(shù)學者都認為筒倉細觀參數(shù)的分布也應該是對稱的，但是實際工程中，筒倉中的顆粒流態(tài)復雜，呈現(xiàn)非均勻性和非對稱性，超壓彈性波的出現(xiàn)和傳播以及倉壁超壓的大小也有微小的差別。以左側倉壁為例，由圖17可見，各監(jiān)測點超壓的出現(xiàn)有先后順序，4號監(jiān)測點的倉壁超壓波峰出現(xiàn)最早，在模擬迭代105 750步時出現(xiàn)穩(wěn)定的成拱超壓，隨著時間推移，4、5、6、7、8 號監(jiān)測點位置的倉壁側壓力波動的波峰依次出現(xiàn)，8號監(jiān)測點超壓波峰出現(xiàn)最晚，在迭代106 510步時才出現(xiàn)穩(wěn)定的成拱超壓。且由圖18可見，雖然右側倉壁彈性波出現(xiàn)的時步與左側倉壁有微小的差別，但在側壁各深度下超壓波動的傳播依然幾乎是均勻的。

6 結論

通過室內實驗及PFC動態(tài)模擬，得出以下結論：

1）卸糧成拱時倉壁出現(xiàn)超壓，導致超壓的原因主要是拱阻礙顆粒出流顆粒速度降低所產生的慣性力。

2）觀察卸糧成拱時倉壁側壓力的波動圖像可以發(fā)現(xiàn)，超壓波峰隨著深度的增加先后出現(xiàn)，說明超壓從拱腳處開始，在倉壁上以波的形式向上傳播，顆粒速度的降低過程也是逐漸向上擴散的。

3）隨著筒倉高度的增加，卸糧成拱時倉壁的側壓力降低，超壓的增量減小，且在筒倉的各深度處超壓系數(shù)也隨深度增加而增大。

4）在筒倉各深度處，出現(xiàn)波峰的時間間隔比較均勻，可以認為超壓波動的傳播也是勻速的。

研究成果揭示了筒倉成拱工況下倉壁超壓的形成機理，相關結論可為筒倉安全性設計提供參考。