基于連續數值模擬的筒倉卸載過程中顆粒物壓強及其速度場分析*

周益嫻

(華北電力大學,非能動核能安全技術北京市重點實驗室,北京 102206)

1 引 言

科學界通常認為顆粒介質(granular media)是大于100 μm的顆粒集合體[1],其布朗效應和顆粒物質之間的凝聚力相較于重力及碰撞力可以忽略.顆粒介質在自然界以及工業生產中廣泛存在,因此對顆粒體系的深入研究將會對全球工業與經濟的發展有極大助益[2].早期顆粒流的研究手段主要是實驗和離散顆粒介質的數值模擬(DEM).離散元數值模擬方法主要是考慮每個個體顆粒物質之間的碰撞[3].該方法發展已較為成熟,但計算速度慢.近年來,法國顆粒物質研究組提出的局部本構理論[4]得到了廣泛使用,這一局部流變理論使得將流體力學方法應用于顆粒流模擬成為可能,這極大程度地降低了運算時間.另外,該理論加深了人們對顆粒流內在物理機理的理解.該方法主要描述的是與固體態相差甚遠的液體狀態,其本質是將顆粒物質看成帶有摩擦性質的非牛頓流體.之后Jop等[5]將此局部流變學從二維向量形式推廣至三維張量形式.該理論被證實能夠應用在許多密集顆粒流構型當中,例如顆粒物坍塌[6]、筒倉卸載[7,8]等.但是存在一定的局限性,第一個局限性是無法描述顆粒物開始流動和停止流動的現象,第二個局限性是無法描述堵塞狀態的運動[9].

顆粒物在筒倉或沙漏中的流動是一個經典的顆粒介質基礎研究現象.人們很早就發現沙漏比水漏作為計時器更加簡便和準確.對于出口在底部的筒倉卸載問題已有大量的研究.當容器尺寸足夠大時,顆粒物的流量始終保持恒定,并且只與出口尺寸有關,這是它的重要特征,但其原因和物理機理仍有待研究.普遍為人們所接受的解釋是Janssen效應[10,11],即容器底部壓強達到一定值后不再隨顆粒高度增加而增加.這是由于顆粒與筒壁間存在著摩擦力,筒壁承擔了部分顆粒重量.由于筒倉出口處的壓強與顆粒物高度無關,因此流量也與高度無關.然而近來一些實驗證明顆粒物的流量與出口處的壓強不具有相關性[12,13].因此Janssen效應并不能完整地解釋此現象.流量Q隨孔徑D變化的經驗公式,早在60年代就由Beverloo 等[14]利用量綱分析和實驗方法給出,并被大量文獻證實相當精確可靠.近幾年人們主要開展了筒倉傾斜情況下的研究[15?17],即出口方向與重力不平行,研究發現流量也保持恒定,并得到了與傾斜角度相關的流量經驗公式.此外Zhou等[18]研究了出口位置在側面(即與地面垂直,傾角為90°)的顆粒物卸載情況(實驗所用筒倉如圖1所示),該工作重點研究了容器厚度W和出口高度D對顆粒物流量的影響.該工作表明,出口在側面的情況下,顆粒物流量仍舊保持恒定,并且作者發現了與出口在底部情況不同的流量規律.當D/W足夠小時,流量規律與出口在底部時的情況一樣,滿足Beverloo經驗公式,即流量為出口面積; 而當足夠大時,流量規律則為可見,在容器厚度很小時,卸載流量呈現了與出口在底部情況下不同的規律.作者使用上述基于局部本構模型的連續模擬方法得到了與實驗結果相同的流量規律,并且證明了導致該結果的主要原因是出口處流線的傾斜角受前后容器壁的摩擦力影響很大.可見,上述連續模擬方法能夠較好地模擬筒倉卸載實驗.本文將在已有工作的基礎上,利用上述基于本構模型的連續數值模擬方法深入分析筒倉卸載現象.首先,分析距離出口近處壓強和速度與W及D之間的關系.進而,由此研究壓強與速度之間的關系,并利用上述數值模擬解釋為何出口在底部時流量規律始終為

2 連續數值模擬

本文采用基于局部本構模型的連續數值模擬方法,即將顆粒物看成非牛頓流體,并且考慮顆粒物之間的摩擦力對顆粒流動的影響.非牛頓不可壓縮的納維-斯托克斯公式可以寫為:

其中D是應變速率張量根據文獻[5]中提出的理論,η為有效黏性系數,它與摩擦系數μ(I)以及局部壓強p相關.D2為應變速率張量的第二不變量,

取流變學常數μs=0.4 ,?μ=0.28和I0=0.4 ,但不考慮文獻[5]提出的體積分數?與I的關系式,此處取?=1.方程(2)中fw代表前后兩個壁面的總體摩擦力.本文實際采用的是二維模擬,為考慮容器前后壁對顆粒物的摩擦以及模擬三維效果,像赫爾-肖(Hele-Shaw)流動一樣,將動量方程沿著y方向進行平均[19,20].根據顆粒流的特性,假設速度沿著y方向幾乎保持不變,因此動量方程中的非線性項前的系數為 1.而黏滯力沿著y方向的積分之和則等于前后壁對顆粒物的摩擦力,假設該力為庫侖力(在每個壁上值為?μwp).該力的方向應當與顆粒物運動方向相反,因此前后兩個壁面的總體摩擦力大小為

取μw=0.1 ,通過改變容器寬度W的值來改變摩擦力的大小.可見本文中容器厚度W的影響等效于容器前后壁摩擦力大小的影響.本文利用上述類三維方法模擬圖1所示實驗,具體計算方法參見文獻[18]和文獻[21].在固體邊界加上非滑移(noslip)的邊界條件,在出口處壓強設為零.數值模擬中容器及顆粒物的尺寸如表1所示(容器相關尺寸參數的定義見圖1).

表1 容器及顆粒物的尺寸Table 1.Size of silo and particle.

3 模擬結果分析

3.1 出口在底部

可見在靠近出口處,動能項開始逐漸增大并且可以部分補償壓強降低,但由于摩擦力的存在,二者之間的關系并不滿足伯努利方程.接下來重點分析距離出口較近處的顆粒物壓強pp和速度vp與容器厚度W及出口高度D的關系.

3.1.1 壓強結果分析

由前面分析可知,顆粒物在距離出口近處壓強和速度變化明顯,在遠離出口的區域內壓強變化較小,速度基本保持恒定,因此本文主要對距離出口近的區域進行分析.下面首先取該區域一個微元薄層進行受力分析.定義一個高度為?λ的微元薄層區域(見圖3中的藍色虛線區域),且若取容器后壁面出口中點處為坐標原點,該區域坐標滿足x∈[?D/2,D/2],y∈[0,W] ,z∈[λ,λ+?λ].假設這個區域的邊界處所受摩擦力與重力平衡,可得:

圖3 出口在底部情況下容器內不同區域受力圖Fig.3.Force diagram of different zones within the silo for the case with orifice at the bottom.

在y和z方向上積分,得到為

在x和z方向積分,可得為

由于重量P=ρgD?λW,取常數α1和α2,方程(5)可寫成

對顆粒物壓強進行無量綱化操作，可得到：

其中γ1=1/(2α1μs),γ2=(α2μw)/(α1μs).(9)式表明顆粒物在靠近出口處壓強僅僅與D/W相關.我們認為z=D處且位于容器中心線上的點也屬于該區域,圖4(a)顯示了不同W值時z=D處無量綱化壓強pp(D)/(ρgL)隨D/L的變化曲線(模擬過程中容器寬度L始終不變).可見,當其他物理量不變時,D增大,出口處壓強隨之增大,這與方程(9)的增減性相符合.且隨著W降低,壓強與D的相關性逐漸降低.當其他物理量不變,W變化時,顆粒物壓強與W強相關.顯見當W降低,顆粒物在z=D處的壓強降低,同樣與方程(9)的增減性相符合.圖4(b)顯示無量綱化壓強pp(D)/(ρgD)隨D/W的變化曲線,模擬數據幾乎重疊并落到一條曲線上,圖中黑色虛線是利用最小二乘法將模擬結果與方程(9)擬合得到的,可見模擬結果與方程(9)符合良好.可得到如下關系:

由方程(10)中γ2=0.49 可推出,當D/W?2時,可得到pp(D)∝ρgD,壓強僅與D相關; 而當D/W?2時,可得到pp(D)∝ρgW,壓強只與W相關.因此根據D/W的大小,壓強呈現出兩種截然不同的相關性.可見出口在底部情況下出口附近的壓強與Zhou 等[18]得到的出口在側面的顆粒物流量有相似的規律.

圖4 在 W 不同的情況下距離出口 D 處的顆粒物壓強pp(D)結果(a)無量綱化壓強 pp(D)/(ρgL)隨無量綱化出口尺寸 D/L 的變化;(b)無量綱化壓強 pp(D)/(ρgD)隨無量綱化出口尺寸 D/W 的變化Fig.4.Results of granular pressure at the distance D from the orifice pp(D)for various W :(a)Dimensionless pressure pp(D)/(ρgL)vs dimensionless orifice size D/L ;(b)dimensionless pressure pp(D)/(ρgD)vs dimensionless orifice size D/W.

3.1.2 出口速度結果分析

大量的實驗表明[22,23],對于出口在底部的情況,流量始終滿足Beverloo方程,可見出口處速度沒有像壓強一樣與相關.文獻[22,23]指出,顆粒物流量與出口中心處法向速度呈正相關,因此圖5(a)給出了無量綱出口中心處法向速度(即豎直方向速度)隨著D/L的變化曲線.由圖5(a)可見,當D/L值較大時,速度隨著W變化會發生輕微改變.圖5(b)畫出了無量綱化出口中心處豎直速度隨無量綱化出口尺寸D/W的變化,圖中黑色虛線為通過最小二乘法得到的方程,可推出出口中心處法向速度與容器厚度及出口尺寸滿足如下關系:

由分母上D/W前√的系數為0.02可知,當D/W?50時,可得到此時流量只與W相關; 而當D/W?50 時,可得此時流量只與D相關,滿足Beverloo定律.實際情況下,D/W?50很難被滿足,因此,在文獻[22,23]中,出口在底部的矩形筒倉卸載問題仍舊滿足Beverloo定理,沒有出現流量與W相關的情況.

圖5 在 W 不同的情況下位于出口處中心線的豎直方向速度 結果(a)無量綱化速度/(2gL)隨無量綱化出口尺寸 D/L 的變化;(b)無量綱化速度/(2gW)隨無量綱化出口尺寸 D/W 的變化Fig.5.Results of vertical velocity on the central streamline for various W :(a)Dimensionless vertical velocity /(2gL)vs dimensionless orifice size D/L ;(b)dimensionless vertical velocity /(2gW)vs dimensionless orifice size D/W.

3.2 出口在側面

圖6顯示了出口在側面的情況下流量恒定期間某一時刻壓強、速度以及無量綱常數I的分布圖.與出口在底部的情況非常相似,壓強在出口處降低而速度顯著增大,I的值在出口處最大,且最大值為 0.1 左右.接下來分析不同區域內的壓強pp和速度vp值與容器厚度及出口尺寸的關系.

3.2.1 壓強結果分析

如圖7所示,研究沿著流線方向距離出口約為D的區域,可假設顆粒物幾乎沿著豎直方向運動,且流動區域沿x方向的尺寸約為D,同出口在底部情況做相同的受力分析,可以得到沿中心流線方向距離出口D處的壓強與D和W的關系為

其中γ1=1/(2α1μs),γ2=(α2μw)/(α1μs).圖8(a)給出了在W不同的情況下距離出口D處的顆粒物無量綱化壓強pp(D)/(ρgL)隨無量綱化出口尺寸D/L的變化,與出口在底部時一樣,當W降低,顆粒物在z=D處的壓強顯著降低.隨著W降低,壓強與D的相關性則逐漸減弱.圖8(b)顯示了無量綱化壓強pp(D)/(ρgD)隨無量綱化出口尺寸D/W的變化,圖中黑色虛線是利用最小二乘法將模擬數據與方程(12)擬合得到的,所得系數與出口在底部的系數非常接近(見圖4(b)),可見,出口處的壓強大小對于出口的位置并不敏感.

3.2.2 出口速度結果分析

對于出口在側面的情況,顆粒物流動方向會發生偏離,將出口處運動方向與豎直方向的夾角記為θ.圖9(a)為無量綱化總速度隨無量綱化出口尺寸D/L的變化,可以看出,出口處總速度大小幾乎與W無關(u0表示出口處水平方向的速度,v0表示豎直方向的速度),且水平速度可被表示為一個與總速度U0以及出口處運動傾斜角相關的函數(θ=arctan(u0/v0)):

圖9(b)顯示了速度傾斜角 sinθ隨無量綱化出口尺寸D/W的變化,圖中黑色虛線為將模擬數據擬合如下方程所得到,可見 sin(θ)與W呈現如下相關性:

圖7 出口在側面情況下容器內不同區域受力圖,圖中陰影部分代表顆粒物停滯區域Fig.7.Force diagram of different zones within the silo for the case with lateral orifice,the dashed area represents the stagnant zone.

上述結果表明,出口處總速度大小只與D有關,但出口速度傾斜角θ與W相關,使得出口中心處法向速度(即水平速度u0)受到W影響,出現了新的流量規律.該現象可從N-S方程中得到解釋,當W減小,摩擦力增大時,速度方向會與重力加速度方向更趨于一致[18].而由方程(14)中分母上D/W前的系數可知,當D/W?4時,可得到 sin(θ)∝W/D,則出口中心處法向速度滿足當D/W?4時,可得到 sin(θ)為常數,則出口中心處法向速度滿足可以推測,對于出口在底部的情況,出口處速度方向基本只與D相關,除非D/W?50.由此可見,側面開口情況下,出口中心處法向速度與W和D的兩種相關性階段的臨界D/W值要遠遠小于出口在底部的情況.

圖8 在 W 不同的情況下距離出口 D 處的顆粒物壓強pp(D)結果(a)無量綱化壓強 pp(D)/(ρgL)隨無量綱化出口尺寸 D/L 的變化;(b)無量綱化壓強 pp(D)/(ρgD)隨無量綱化出口尺寸 D/W 的變化Fig.8.Results of granular pressure at the distance D from the orifice pp(D)for various W :(a)Dimensionless pressure pp(D)/(ρgL)vs.dimensionless orifice size D/L ;(b)dimensionless pressure pp(D)/(ρgD)vs.dimensionless orifice size D/W.

圖9 在 W 不同的情況下位于出口處中心流線上速度結果(a)無量綱化總速度(2gL)隨無量綱化出口尺寸D/L 的變化;(b)速度傾斜角 sinθ 隨無量綱化出口尺寸D/W的變化Fig.9.Results of velocity on the central streamline at the orifice for various W :(a)Dimensionless total velocity /(2gL)vs dimensionless orifice size D/L ;(b)angle of inclination sinθ vs dimensionless orifice size D/W.

4 結 論

本文主要通過基于局部本構理論的連續數值模擬方法研究出口在底部和側面的筒倉卸載問題.深入研究了容器厚度W(影響容器前后壁摩擦力大小)和出口尺寸D對顆粒物內部物理量的影響,著重分析了靠近出口處的壓強及速度.首先通過受力分析預測了不同位置的壓強與W以及D的關系,并利用模擬方法進行了驗證.理論分析和數值模擬結果都表明,距離出口較近區域的壓強存在兩種相關性,當D/W足夠小時,壓強只與D相關,當D/W足夠大時,壓強只與W相關.出口在底部和側面兩種情況均呈現上述規律.隨后,分析了出口處法向速度與W以及D的相關性,有趣的是,出口中心處法向速度呈現出與壓強不同的相關性.對于出口在底部的情況,對于模擬中任意D/W值,出口中心處法向速度仍然只與D相關; 而出口在側面時,當D/W足夠小時,出口中心處法向速度只與D相關,當D/W足夠大時,出口中心處法向速度只與W相關.這個結果與一直以來人們所認為的壓強是影響出口速度的主要因素明顯相違背[10],而與近來Perge等[12]以及Aguirre等[13]所做實驗得出的結論一致,即筒倉出口處顆粒物流量與壓強無關.對于出口在側面的情況,總速度大小只與D相關,而出口中心處法向速度(水平速度)與W相關,其原因是出口處速度傾斜角與W相關,當W降低,出口處流線會更靠近豎直方向而總速度大小不變,因此出口中心處法向速度(即水平方向速度)降低.對于出口在底部的矩形筒倉卸載問題,本文還揭示了當容器厚度W與出口高度D相比很小時,流量仍然滿足Beverloo定律的原因[22,23].這是由于與出口在側面相比,出口在底部時,造成相關性規律改變的D/W臨界值較大,只有當D/W?50 時,出口流速(豎直方向速度)才會受W影響,呈現不同的變化規律,一般實際情況無法滿足該條件.

綜上所述,本文利用基于局部本構理論的連續數值模擬得到了與文獻中實驗一致的結果[18,22,23],并深入剖析了顆粒物內部運動機制,可見該模擬方法能夠較好地預測顆粒流運動規律.從前文分析中可看出,停滯區的角度對流量影響較大,而Kamrin和Koval[24]所提出的非局部理論能夠更好地描述該區域的運動情況,未來可將局部本構模型修正為更為精確的非局部本構模型,并比較二者在筒倉卸載問題中結果的區別.