基于改進(jìn)顆粒模型的筒倉卸糧宏細(xì)觀力學(xué)響應(yīng)模擬研究

馮 永 潘 樊 劉 杰

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

筒倉在貯存散體物料方面具有容量大、占地面積小等優(yōu)點(diǎn)而被廣泛應(yīng)用，其結(jié)構(gòu)設(shè)計(jì)的關(guān)鍵在于對倉儲工況荷載的正確把握[1, 2]。而相對于靜態(tài)工況，糧食等散體物料在筒倉內(nèi)的動態(tài)運(yùn)動規(guī)律則更為復(fù)雜[3-8]，較多學(xué)者對此展開研究。

圖1 顆粒接觸模型

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

本研究在單一ball單元離散元程序基礎(chǔ)上，添加黏度系數(shù)較大的微型顆粒來模擬粉塵，旨在建立一種改進(jìn)顆粒組構(gòu)力學(xué)模型，并將該模型應(yīng)用于筒倉卸糧模擬研究中?；谛屡f模型的對比，對卸糧過程中散體動力學(xué)參數(shù)進(jìn)行深入探究和分析，以期為改進(jìn)模型、筒倉結(jié)構(gòu)優(yōu)化提供參考。

1 模型本構(gòu)關(guān)系

1.1 傳統(tǒng)ball單元模型與改進(jìn)顆粒模型比較分析

根據(jù)離散元法[19]，傳統(tǒng)ball單元模型采用單一圓形顆粒模擬糧食顆粒，顆粒單元用i表示，顆粒組構(gòu)見圖2a?；谀壳皞鹘y(tǒng)ball單元顆粒模型的局限性，本研究的改進(jìn)顆粒模型將傳統(tǒng)ball單元顆粒改成2種大小不同的顆粒，新加入的顆粒單元模擬雜質(zhì)用k表示，改進(jìn)顆粒模型包括2種基本單元大球顆粒i和小球顆粒k，改進(jìn)顆粒組構(gòu)如圖2b所示。

注：Kn為顆粒正向剛度/N/m；Ks為顆粒切向剛度/N/m；Tn為正向強(qiáng)度/N/m2；Ts為切向強(qiáng)度/N/m2；u為摩擦系數(shù)；i為大型顆粒單元體，k為小型顆粒單元體。圖2 傳統(tǒng)ball單元模型與改進(jìn)顆粒模型

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

1.2 改進(jìn)顆粒模型基本單元間接觸關(guān)系

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

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

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

(1)

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

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

2 物理模型實(shí)驗(yàn)及離散元模型的建立

2.1 實(shí)驗(yàn)裝置

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

將實(shí)驗(yàn)?zāi)Ｐ蛡}內(nèi)裝滿粒徑為5～6 mm的小麥，通過室內(nèi)實(shí)驗(yàn)測得實(shí)驗(yàn)小麥的具體參數(shù)，見表1，待裝料穩(wěn)定后獲取筒倉靜態(tài)側(cè)壓力，之后撤掉料斗底板進(jìn)行卸糧實(shí)驗(yàn)獲取動態(tài)側(cè)壓力，多次實(shí)驗(yàn)取其平均值。

表1 實(shí)驗(yàn)小麥參數(shù)

2.2 實(shí)驗(yàn)結(jié)果

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

圖4 實(shí)驗(yàn)裝置

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

圖5 倉壁側(cè)壓力隨時間變化圖

2.3 離散元模型及參數(shù)的建立

為了深入探究改進(jìn)模型在筒倉模擬研究中的應(yīng)用優(yōu)勢，參照室內(nèi)實(shí)驗(yàn)筒倉，建立同尺寸筒倉模型。

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

表2 傳統(tǒng)ball單元顆粒及改進(jìn)后顆粒參數(shù)

表3 計(jì)算時步參數(shù)

2.4 驗(yàn)證模擬模型的合理性

筒倉靜態(tài)側(cè)壓力的計(jì)算采用目前普遍應(yīng)用的Janssen理論，計(jì)算公式為：

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

(2)

k=1.1(1-sinφ)

(3)

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

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

采用Z-Score標(biāo)準(zhǔn)化對實(shí)驗(yàn)數(shù)據(jù)和2次模擬數(shù)據(jù)進(jìn)行了無量綱處理，如圖6c，模擬值與實(shí)驗(yàn)值整體契合度較好，其中3號點(diǎn)位處的偏差最大，傳統(tǒng)模擬和改進(jìn)模擬分別相對實(shí)驗(yàn)值偏離0.20、0.39個標(biāo)準(zhǔn)差，其余點(diǎn)位的偏離值均在0.18以內(nèi)，與實(shí)驗(yàn)契合度高，可以驗(yàn)證模擬的合理性。

圖6 倉壁靜態(tài)側(cè)壓力對比圖

3 倉壁壓力數(shù)據(jù)分析

3.1 壓力曲線趨勢對比分析

筒倉在自由卸糧過程中，位于上部的10號監(jiān)測點(diǎn)測得的動態(tài)側(cè)壓力曲線波動幅度較小，隨著糧食的下卸迅速進(jìn)入零壓力階段，故選取2、4、6、8這4個監(jiān)測點(diǎn)的卸糧動態(tài)側(cè)壓力來進(jìn)行數(shù)據(jù)分析。

圖7為筒倉各監(jiān)測點(diǎn)的動態(tài)側(cè)壓力模擬對比圖，在卸糧過程中，倉壁各個監(jiān)測點(diǎn)的動態(tài)側(cè)壓力均大于靜止?fàn)顟B(tài)下的側(cè)壓力；隨著顆粒的不斷流出，各監(jiān)測點(diǎn)的動態(tài)側(cè)壓力上下波動，且隨著時步的增加，側(cè)壓力逐漸減小至零。

圖7 各點(diǎn)位動態(tài)側(cè)壓力模擬對比圖

對比圖8中2種模擬結(jié)果可以發(fā)現(xiàn)，顆粒單元改進(jìn)后卸料所需時步遠(yuǎn)遠(yuǎn)大于改進(jìn)前。顆粒單元改進(jìn)前的倉壁動態(tài)側(cè)壓力曲線波動幅度較小，各監(jiān)測點(diǎn)處的壓力曲線交叉較少；而改進(jìn)后的模擬曲線波動幅度較大，各監(jiān)測點(diǎn)處的側(cè)壓力曲線交叉較多。

圖8 改進(jìn)模型與傳統(tǒng)模型倉壁動態(tài)側(cè)壓力整體對比圖

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

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

3.2 數(shù)據(jù)準(zhǔn)確性對比分析

糧食對倉壁的動態(tài)側(cè)壓力，可根據(jù)GB 50077—2017《鋼筋混凝土筒倉設(shè)計(jì)標(biāo)準(zhǔn)》[22]進(jìn)行計(jì)算。

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

(4)

k=tan2(45°-φ/2)

(5)

式中：Ch為深倉貯料水平壓力修正系數(shù)，根據(jù)GB 50077—2017《鋼筋混凝土筒倉設(shè)計(jì)標(biāo)準(zhǔn)》計(jì)算:監(jiān)測點(diǎn)距倉壁底部的深度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 不同方法下測得的動態(tài)側(cè)壓力

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

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

(6)

圖10 模擬與實(shí)驗(yàn)過程中監(jiān)測點(diǎn)切應(yīng)力與位移關(guān)系曲線圖

4 倉壁壓力細(xì)觀力學(xué)機(jī)理分析

在細(xì)觀角度上，通過卸糧過程顆粒的力與位移的變化關(guān)系以及力鏈的構(gòu)成，深入探究倉壁動態(tài)側(cè)壓力變化的機(jī)理。取改進(jìn)前、后顆粒模型的2、4、6這3個監(jiān)測點(diǎn)，做切應(yīng)力與位移曲線對比圖，如圖10所示。傳統(tǒng)ball單元模型在卸糧的過程中，切應(yīng)力的增長相對緩慢，說明顆粒間黏結(jié)強(qiáng)度較低。而改進(jìn)顆粒模型各監(jiān)測點(diǎn)切應(yīng)力迅速提高，并很快達(dá)到最大值，說明顆粒間黏結(jié)強(qiáng)度較高。取實(shí)驗(yàn)與模擬的6號監(jiān)測點(diǎn)剪切力-位移曲線，以實(shí)驗(yàn)?zāi)Ｐ颓€為標(biāo)準(zhǔn)值，根據(jù)式(6)，得出傳統(tǒng)ball單元模型曲線的標(biāo)準(zhǔn)差為100.1 kPa，改進(jìn)顆粒模型曲線的標(biāo)準(zhǔn)差為87.5 kPa，比傳統(tǒng)ball單元曲線標(biāo)準(zhǔn)差減小12.5%。改進(jìn)顆粒模型所得數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)更接近。

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

圖11 接觸力鏈圖

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

從細(xì)觀角度來看，傳統(tǒng)ball單元模型由于切應(yīng)力提高較慢，顆粒內(nèi)部結(jié)構(gòu)應(yīng)力在不斷地調(diào)整，顆粒間摩擦力、黏聚力較小，卸糧過程比較順暢，無法形成拱效應(yīng)。而改進(jìn)顆粒模型由于切應(yīng)力迅速提高，顆粒之間的裂隙被壓實(shí)，剪切帶完全貫通，導(dǎo)致顆粒之間的黏聚力增加，摩擦更為充分，卸糧過程中出現(xiàn)拱效應(yīng)，更符合實(shí)驗(yàn)過程。

5 結(jié)論

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

以實(shí)驗(yàn)倉壁側(cè)壓力曲線為標(biāo)準(zhǔn)值，傳統(tǒng)ball單元模型的側(cè)壓力曲線的標(biāo)準(zhǔn)差為606.66 Pa，改進(jìn)顆粒模型曲線的標(biāo)準(zhǔn)差為137.50 Pa，相比減小77.3%，改進(jìn)顆粒模型數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)更接近。

以實(shí)驗(yàn)剪切力-位移曲線為標(biāo)準(zhǔn)值，傳統(tǒng)ball單元模型剪切力與位移曲線的標(biāo)準(zhǔn)差為100 kPa，改進(jìn)顆粒模型的標(biāo)準(zhǔn)差為87.5 kPa，相比減小12.5%，擬合的效果更好。