基于改進顆粒模型的筒倉卸糧宏細觀力學響應模擬研究

馮 永 潘 樊 劉 杰

(河南工業大學土木工程學院，鄭州 450001)

筒倉在貯存散體物料方面具有容量大、占地面積小等優點而被廣泛應用，其結構設計的關鍵在于對倉儲工況荷載的正確把握[1, 2]。而相對于靜態工況，糧食等散體物料在筒倉內的動態運動規律則更為復雜[3-8]，較多學者對此展開研究。

圖1 顆粒接觸模型

陳長冰等[9]用PFC2D中的等圓形單元模擬顆粒，研究了筒倉卸料過程中儲料的流動狀態、顆粒間接觸力以及側壁壓力的變化。丁盛威等[10]運用PFC3D中的球形單元，分別對靜態工況及動態工況的筒倉倉壁側壓力進行模擬分析，得出顆粒流動狀態影響倉壁壓力變化。有研究通過離散元模擬和室內實驗分析了倉壁內部附著物對側壓力分布的影響，以及動態側壓力分布規律與物料流動狀態，但其模擬所采用的顆粒接觸仍是單一球單元的線性接觸[11-14]。隨著對側壓力研究的不斷深入，不少學者提出筒倉成拱效應，有研究通過結拱時候的拱線分析,得出糧食速度驟降而產生慣性力進一步增大了倉壁側壓,并通過實驗及數值模擬進行了驗證分析，但所采用的的模擬單元依然是傳統ball單元[4, 15-17]。已有報道雖然一定程度上推動了筒倉內顆粒流動模擬的研究，但其采用的模擬單元模型并不十分嚴謹，實際中在糧食周圍存在的雜質及粉塵對顆粒流動及側壓力均有很大影響。現有研究采用的單一ball單元模型，簡圖如圖1所示，由于顆粒粒徑相同，顆粒在排列時空隙較大、擠壓不充分，不能客觀反映顆粒間及顆粒與倉壁間接觸的實際情況，特別是顆粒-倉壁接觸面積變小[18]，常導致數值模擬結果和實驗結果不吻合。

本研究在單一ball單元離散元程序基礎上，添加黏度系數較大的微型顆粒來模擬粉塵，旨在建立一種改進顆粒組構力學模型，并將該模型應用于筒倉卸糧模擬研究中。基于新舊模型的對比，對卸糧過程中散體動力學參數進行深入探究和分析，以期為改進模型、筒倉結構優化提供參考。

1 模型本構關系

1.1 傳統ball單元模型與改進顆粒模型比較分析

根據離散元法[19]，傳統ball單元模型采用單一圓形顆粒模擬糧食顆粒，顆粒單元用i表示，顆粒組構見圖2a。基于目前傳統ball單元顆粒模型的局限性，本研究的改進顆粒模型將傳統ball單元顆粒改成2種大小不同的顆粒，新加入的顆粒單元模擬雜質用k表示，改進顆粒模型包括2種基本單元大球顆粒i和小球顆粒k，改進顆粒組構如圖2b所示。

注：Kn為顆粒正向剛度/N/m；Ks為顆粒切向剛度/N/m；Tn為正向強度/N/m2；Ts為切向強度/N/m2；u為摩擦系數；i為大型顆粒單元體，k為小型顆粒單元體。圖2 傳統ball單元模型與改進顆粒模型

圖2中的彈簧單元表示顆粒間接觸前的線彈性行為，阻尼器單元用來產生能量耗散。顆粒間未接觸時無作用力，用分離器單元來表示。顆粒切向接觸的滑片單元遵循摩爾-庫倫摩擦定律，切應力大小應由顆粒間的黏度系數與摩擦系數共同確定。

1.2 改進顆粒模型基本單元間接觸關系

改進顆粒模型單元間的接觸，包括顆粒-顆粒間的接觸與顆粒-倉壁間的接觸，2種接觸情況計算可選用Lin等[20]提出的幾何平均法，解出顆粒間或顆粒與倉壁間的接觸點坐標，并將其轉化成剪切力與位移關系曲線，顆粒與倉壁間接觸關系如圖3所示。

注：xci，yci；xck，yck為圓心坐標/mm；xc，yc為接觸點坐標/mm；xcj，ycj為倉壁j起始坐標/mm；ri，rk為接觸半徑的向量/mm；θi，θk為角速度/rad/s；vxi、vyi，vxk、vyk為線速度為交點向量，αi、αk為與X軸夾角/rad，rw為接觸半徑的向量。圖3 顆粒與倉壁接觸示意圖

假設坐標系中半徑為ri的圓形顆粒i的圓心坐標為(xci，yci)，倉壁j的起始坐標為(xcj,ycj)，橫軸與x軸的正向夾角為θi。圓形顆粒i的方程為：

(1)

若式(1)無解，則顆粒與倉壁間無接觸；若有解，則證明顆粒與倉壁有接觸，此解即為接觸點坐標(xc，yc)。

基本的單元間力與位移的關系在已有工作中已經證明[21]，本研究不再贅述。

2 物理模型實驗及離散元模型的建立

2.1 實驗裝置

為了觀察筒倉內糧食顆粒的流動狀態，采用亞克力材料制作模型筒倉倉壁，筒倉根據實際常用筒倉尺寸1/20的縮尺比例建立，倉壁高1 000 mm,直徑為500 mm，倉壁厚度約為8 mm,漏斗半頂角為60°。按照等間距的原則在倉壁左右兩側布置監測點來檢測筒倉卸糧過程中的側壓力。筒倉模型的具體尺寸和監測點布置如圖4所示。

將實驗模型倉內裝滿粒徑為5～6 mm的小麥，通過室內實驗測得實驗小麥的具體參數，見表1，待裝料穩定后獲取筒倉靜態側壓力，之后撤掉料斗底板進行卸糧實驗獲取動態側壓力，多次實驗取其平均值。

表1 實驗小麥參數

2.2 實驗結果

筒倉卸糧實驗進行中，由于顆粒初始流動速度較快，最上部監測點的數據差異性較大，故在數據提取處理時選取2、4、6、8，這4個監測點位的數據，通過采集的倉壁應力數據，整理得出倉壁側壓力隨時間變化圖，如圖4所示。

圖4 實驗裝置

由圖5可以看出，隨著卸料的進行，4個點位的側壓力曲線走向基本保持一致，倉壁壓力隨實踐步長的增加逐漸趨向于零。倉壁的側壓力呈現不規律性的振蕩，且越接近筒倉底部，倉壁側壓力振蕩幅度越大。

圖5 倉壁側壓力隨時間變化圖

2.3 離散元模型及參數的建立

為了深入探究改進模型在筒倉模擬研究中的應用優勢，參照室內實驗筒倉，建立同尺寸筒倉模型。

改進顆粒模型中加入體積更小但黏度更大的顆粒k，使i、k顆粒數量比為10∶1，用來模擬筒倉中糧食顆粒間的粉塵雜質。模擬的顆粒及模型時步參數如表2、表3所示。

表2 傳統ball單元顆粒及改進后顆粒參數

表3 計算時步參數

2.4 驗證模擬模型的合理性

筒倉靜態側壓力的計算采用目前普遍應用的Janssen理論，計算公式為：

Ph=γρ(1-e-μks/ρ)/μ

(2)

k=1.1(1-sinφ)

(3)

式中：γ為儲料的重力密度/kN/m3；ρ為筒倉水平凈截面的水力半徑/m；μ為儲料與倉壁的摩擦系數；k為側壓力系數；s為儲料頂面至所計算截面的距離/m；φ為儲料的內摩擦角/°。

選取改進前、后模型倉壁右側的靜態儲糧側壓力數據、物理實驗的靜態儲糧側壓力數據、式(2)～式(3)計算的靜態側壓力值繪成曲線圖(圖6)。結果顯示，Janssen公式計算所得的靜態側壓力曲線與實驗狀態下曲線基本吻合，圖6b為兩者的誤差線，其中最大相對誤差出現在1號點位處，為14.6%，7號點位的相對誤差為10.1%，其余點位誤差都低于9%。

采用Z-Score標準化對實驗數據和2次模擬數據進行了無量綱處理，如圖6c，模擬值與實驗值整體契合度較好，其中3號點位處的偏差最大，傳統模擬和改進模擬分別相對實驗值偏離0.20、0.39個標準差，其余點位的偏離值均在0.18以內，與實驗契合度高，可以驗證模擬的合理性。

圖6 倉壁靜態側壓力對比圖

3 倉壁壓力數據分析

3.1 壓力曲線趨勢對比分析

筒倉在自由卸糧過程中，位于上部的10號監測點測得的動態側壓力曲線波動幅度較小，隨著糧食的下卸迅速進入零壓力階段，故選取2、4、6、8這4個監測點的卸糧動態側壓力來進行數據分析。

圖7為筒倉各監測點的動態側壓力模擬對比圖，在卸糧過程中，倉壁各個監測點的動態側壓力均大于靜止狀態下的側壓力；隨著顆粒的不斷流出，各監測點的動態側壓力上下波動，且隨著時步的增加，側壓力逐漸減小至零。

圖7 各點位動態側壓力模擬對比圖

對比圖8中2種模擬結果可以發現，顆粒單元改進后卸料所需時步遠遠大于改進前。顆粒單元改進前的倉壁動態側壓力曲線波動幅度較小，各監測點處的壓力曲線交叉較少；而改進后的模擬曲線波動幅度較大，各監測點處的側壓力曲線交叉較多。

圖8 改進模型與傳統模型倉壁動態側壓力整體對比圖

其中，改進前，2號和8號監測點側壓力降為零的時步分別為252、17.5萬步；2號監測點的最大側壓力在第105萬步時為3 177 Pa，8號監測點的最大側壓力在第63萬步時為1 415 Pa。改進后，2號和8號監測點處的側壓力減小為零的時步增加至644、420萬步；2號監測點最大側壓力在第15.4萬步時為4 364 Pa ；8號監測點處最大側壓力在第84萬步時為2 292 Pa。相比改進前模型，2號和8號監測點時步數分別增加了155%和140%，終止時步數明顯增大；最大側壓力分別增大了37.4%和62%，振蕩的數值范圍也變大。

主要原因是改進前的顆粒單元是單一球形顆粒，顆粒間空隙較大，接觸擠壓不充分又缺少黏結力，所以卸料過程比較順暢，沒有出現明顯的壓力成拱現象，即倉壁動態壓力波動較小，各監測點壓力曲線交叉較少且基本呈現出整體流動狀態，直到顆粒完全流出，這一過程卸糧時間較短；改進后的模型添加了黏度大、直徑小的圓形顆粒模擬粉塵，筒倉內顆粒分布更緊密，顆粒間的黏結、摩擦力使得顆粒下卸時出現壓力拱，在倉壁處壓力急速增大，隨著顆粒的不斷流動，壓力拱被破壞，倉壁壓力又迅速變小，周而復始的成拱-破拱使得倉壁側壓力上下波動且幅度較大。又因為壓力拱在卸糧過程位置不斷的變化，所以各個監測點處壓力的數值都會很大，各點位的壓力曲線交叉較多。同樣基于此，倉壁兩側的顆粒流動受到阻塞，顆粒完全卸出用時較長。

3.2 數據準確性對比分析

糧食對倉壁的動態側壓力，可根據GB 50077—2017《鋼筋混凝土筒倉設計標準》[22]進行計算。

Ph=Chγρ(1-e-μks/ρ/μ)

(4)

k=tan2(45°-φ/2)

(5)

式中：Ch為深倉貯料水平壓力修正系數，根據GB 50077—2017《鋼筋混凝土筒倉設計標準》計算:監測點距倉壁底部的深度h為0～0.66 m時，Ch=2；h=0.72時，Ch=1.82；h=0.81時，Ch=1.56；h=0.9，Ch=1.29。

圖9 不同方法下測得的動態側壓力

圖9為運用不同方法下的倉壁動態側壓力對比圖。由規范計算公式得到的各測點動態側壓力隨距倉底高度的增加呈現出平穩減小的趨勢，且在下降過程中無波動。但實驗數據和改進后的模擬曲線顯示，在2號和4號點位處動態側壓力有向上波動的趨勢，分別超出規范計算數值6.6%和4.1%，表明筒倉在卸糧過程中，在靠近筒倉底部1/3處有成拱效應，糧食在此處堆積導致倉壁動態側壓力較周圍測點有明顯的升高，這與張大英等[23]的研究結果相符。而其他監測點處的動態側壓力顯然比規范計算值更大，這符合規范取值保守的原則，故改進后的模擬數據與實驗數據更相符。

用改進前、后的動態側壓力和實驗數據做對比，改進前的倉壁動態側壓力模擬值與實驗值的相對誤差較大，其絕對值在7.93%～30.06%之間；而改進后的倉壁動態側壓力模擬值與實驗值吻合較好，相對誤差絕對值在1.64%～10.93%之間。改進后各測點模擬曲線與計算曲線相對誤差絕對值均小于11%，在合理范圍，可以參考相近研究的誤差范圍[24]。以圖8的數據曲線為標準值，利用式(6)求出改進前模型曲線的標準差為606.66 Pa，改進后顆粒模型曲線的標準差為137.50 Pa，比改進前模型減小了77.3%，與實驗數據曲線更接近，擬合效果更好。

(6)

圖10 模擬與實驗過程中監測點切應力與位移關系曲線圖

4 倉壁壓力細觀力學機理分析

在細觀角度上，通過卸糧過程顆粒的力與位移的變化關系以及力鏈的構成，深入探究倉壁動態側壓力變化的機理。取改進前、后顆粒模型的2、4、6這3個監測點，做切應力與位移曲線對比圖，如圖10所示。傳統ball單元模型在卸糧的過程中，切應力的增長相對緩慢，說明顆粒間黏結強度較低。而改進顆粒模型各監測點切應力迅速提高，并很快達到最大值，說明顆粒間黏結強度較高。取實驗與模擬的6號監測點剪切力-位移曲線，以實驗模型曲線為標準值，根據式(6)，得出傳統ball單元模型曲線的標準差為100.1 kPa，改進顆粒模型曲線的標準差為87.5 kPa，比傳統ball單元曲線標準差減小12.5%。改進顆粒模型所得數據與實驗數據更接近。

傳統ball單元模型與改進顆粒模型卸糧試樣的接觸力鏈圖如圖11所示，其中右側放大部分較粗的線條為顆粒力的傳遞。

圖11 接觸力鏈圖

由圖11a可以看出，顆粒接觸力方向多為豎直向下,很少出現水平力鏈，這是由于單一顆粒在筒倉卸料過程中顆粒與顆粒的摩擦，以及顆粒與倉壁的摩擦都不充分，無法形成拱效應，因此顆粒對倉壁的擠壓程度輕，水平力鏈較少。在筒倉卸料過程中(圖11b)，顆粒力的傳遞方向包括水平方向和豎直方向，且力鏈分叉現象明顯,這是由于加入的小顆粒k使得大小顆粒受壓后孔隙減小，顆粒排列更為密實，在卸料過程中，顆粒-顆粒、顆粒-倉壁間摩擦充分，在倉壁處出現拱效應，水平力大幅增加。

從細觀角度來看，傳統ball單元模型由于切應力提高較慢，顆粒內部結構應力在不斷地調整，顆粒間摩擦力、黏聚力較小，卸糧過程比較順暢，無法形成拱效應。而改進顆粒模型由于切應力迅速提高，顆粒之間的裂隙被壓實，剪切帶完全貫通，導致顆粒之間的黏聚力增加，摩擦更為充分，卸糧過程中出現拱效應，更符合實驗過程。

5 結論

針對傳統ball單元顆粒模型對筒倉卸糧模擬研究的局限性，本研究提出了改進顆粒本構模型的思想，并運用物理數學方法推導了模型的可行性。本研究基于改進模型思想，依托室內實驗，通過離散元軟件進行筒倉卸糧模擬，從宏細觀角度分別對比驗證了改進模擬相較于傳統ball單元模型在模擬中的優勢。結果表明，倉壁側壓力隨著時步的進行逐漸減小至零；2號和10號監測點處相比于傳統ball單元模型，改進顆粒模型的倉壁壓力變至為零的時步增加了155%和140%，完全卸糧時間大幅增加，最大側壓力增大了37.4%和62%，振蕩的數值范圍更大。

以實驗倉壁側壓力曲線為標準值，傳統ball單元模型的側壓力曲線的標準差為606.66 Pa，改進顆粒模型曲線的標準差為137.50 Pa，相比減小77.3%，改進顆粒模型數據與實驗數據更接近。

以實驗剪切力-位移曲線為標準值，傳統ball單元模型剪切力與位移曲線的標準差為100 kPa，改進顆粒模型的標準差為87.5 kPa，相比減小12.5%，擬合的效果更好。