吳 坤 劉向軍 戴椰凌

(北京科技大學(xué)能源與環(huán)境工程學(xué)院,北京 100083)

引言

顆粒移動(dòng)床在工業(yè)生產(chǎn)中具有很廣泛的應(yīng)用,如燃料熱解氣化[1],煙氣除塵[2],高溫散料余熱回收等[3-4].顆粒移動(dòng)床內(nèi)的流動(dòng)特性是工藝設(shè)計(jì)與優(yōu)化運(yùn)行的重要基礎(chǔ),尚待深入地揭示.

現(xiàn)有的描述移動(dòng)床顆粒流的研究方法有離散元方法(DEM)和連續(xù)性假設(shè)方法.離散元方法從單個(gè)顆粒出發(fā)研究每個(gè)顆粒與其他顆粒作用的詳細(xì)經(jīng)歷[5],模擬精度高,近年來(lái)有很多學(xué)者應(yīng)用此方法對(duì)移動(dòng)床內(nèi)顆粒的流動(dòng)情況進(jìn)行研究,包括移動(dòng)床的幾何構(gòu)造和顆粒特性對(duì)顆粒出流量的影響[6-8],移動(dòng)床壁面的應(yīng)力分布[9-10]以及移動(dòng)床內(nèi)顆粒繞管特性[11]等.但離散元方法的計(jì)算量龐大,且對(duì)復(fù)雜形狀的工業(yè)顆粒的表征尚待改進(jìn),對(duì)實(shí)際大型工業(yè)移動(dòng)床的模擬尚存在困難.連續(xù)介質(zhì)假設(shè)方法將移動(dòng)床內(nèi)的顆粒流視為擬流體,在歐拉坐標(biāo)下求解方程.擬流體模型計(jì)算量小,其關(guān)鍵是建立顆粒流的本構(gòu)關(guān)系.

顆粒流一般認(rèn)為有3 種流動(dòng)狀態(tài):顆粒緊密堆積,流動(dòng)非常緩慢,顆粒之間依靠摩擦相互作用的準(zhǔn)靜態(tài)狀態(tài)[12];顆粒非常稀疏,顆粒之間以二元碰撞為主的快速流動(dòng)狀態(tài)[13];以及處于二者之間,摩擦和碰撞對(duì)顆粒流動(dòng)均有影響的中間流動(dòng)狀態(tài)[14-15].這3 種流動(dòng)狀態(tài)也被稱(chēng)為顆粒固態(tài),顆粒氣態(tài)和顆粒液態(tài).顆粒在不同流態(tài)中呈現(xiàn)不同的流動(dòng)特性對(duì)于顆粒擬流體方法的發(fā)展是一大挑戰(zhàn),顆粒固態(tài)時(shí)一般使用黏塑性理論描述,顆粒氣態(tài)時(shí)使用顆粒流動(dòng)理學(xué)理論(KTGF)[16],對(duì)于顆粒液態(tài),還沒(méi)有形成統(tǒng)一的理論去描述,一般可以將準(zhǔn)靜態(tài)的黏塑性理論拓展到相對(duì)稀疏的顆粒流,得到摩擦黏度模型,或者通過(guò)修正顆粒動(dòng)理學(xué)理論實(shí)現(xiàn)對(duì)顆粒液態(tài)的描述[17-18].對(duì)于移動(dòng)床內(nèi)的密集顆粒流,顆粒流存在幾乎不流動(dòng)的滯止區(qū),存在速度梯度很大的剪切區(qū),也存在流動(dòng)較均勻的主流區(qū),不同結(jié)構(gòu)的移動(dòng)床內(nèi)部可能存在2 種或3 種顆粒狀態(tài).現(xiàn)有的普遍認(rèn)為能用于描述這種密集顆粒流動(dòng)的模型有Shaeffer 模型[19],S-S模型[20],μ(I)模型[21]等,其中μ(I)模型基于黏塑性理論,形式簡(jiǎn)單,模型參數(shù)與顆粒物性關(guān)聯(lián)直接,近年來(lái)受到國(guó)內(nèi)外學(xué)者廣泛重視[22].

μ(I)模型將顆粒物料處理成黏塑性流體,引入無(wú)量綱慣性數(shù)I來(lái)描述顆粒的流動(dòng)狀態(tài)與特性,形式簡(jiǎn)單[21].目前μ(I)模型在山體滑坡、泥石流、雪崩等顆粒坍塌的外流流動(dòng)預(yù)測(cè)應(yīng)用廣泛[23-25].近年來(lái)不斷有研究者將μ(I)應(yīng)用于不同結(jié)構(gòu)的移動(dòng)床流動(dòng)模擬,Chauchat 和Medale[26]采用有限元的方法,通過(guò)改進(jìn)算法處理壓力和黏度的強(qiáng)依賴(lài)關(guān)系,并將改進(jìn)后的模型應(yīng)用于斜面顆粒流動(dòng)和顆粒的圓柱繞流,和理論以及實(shí)驗(yàn)的結(jié)果吻合較好;Tian 等[22]對(duì)移動(dòng)床卸料過(guò)程進(jìn)行模擬,并將μ(I)模型與DEM 模擬、Schasffer 模型、S-S 模型進(jìn)行對(duì)比研究;Bartsch等[27]將μ(I)模型應(yīng)用于移動(dòng)床熱交換器,對(duì)顆粒流在排管間的流動(dòng)進(jìn)行了模擬研究;Baumann 等[28]將μ(I)模型與傳熱模擬相結(jié)合研究移動(dòng)床內(nèi)顆粒流與管道的換熱,其重點(diǎn)在于研究傳熱,對(duì)μ(I)模型和流動(dòng)計(jì)算簡(jiǎn)化較大;Schneiderbauer 等[29]將動(dòng)理學(xué)理論模型和μ(I)相結(jié)合,應(yīng)用于移動(dòng)床和流化床顆粒流模擬;Staron 等[30]應(yīng)用μ(I)模型模擬移動(dòng)床卸料過(guò)程并和DEM 方法進(jìn)行對(duì)比.

在上述μ(I)模型的各種應(yīng)用中,一類(lèi)是將顆粒流視為密度、顆粒體積分?jǐn)?shù)不變的不可壓縮單相μ(I)流體,另一類(lèi)則是聯(lián)立氣相控制方程,考慮顆粒體積分?jǐn)?shù)變化的多相流處理方法.對(duì)于移動(dòng)床內(nèi)顆粒流動(dòng)的模擬,采用現(xiàn)有不可壓縮單相μ(I)模型,床內(nèi)顆粒流密度和顆粒體積分?jǐn)?shù)為不變的常數(shù),無(wú)法預(yù)報(bào)顆粒體積分?jǐn)?shù)變化,忽略了流動(dòng)傳熱的關(guān)鍵局部細(xì)節(jié),誤差較大[25],而若采用多相流模型,又將絕大部分的計(jì)算量浪費(fèi)在對(duì)顆粒力學(xué)作用可以忽略不計(jì)的氣相方程求解上.μ(I)模型本質(zhì)上是一種與局部流動(dòng)狀態(tài)密切相關(guān)的流變模型,局部流動(dòng)與局部壓力、顆粒體積分?jǐn)?shù)密切相關(guān)[31].本文基于單相μ(I)模型,補(bǔ)充局部顆粒體積分?jǐn)?shù)與顆粒局部壓力和局部似密度的關(guān)系式,將移動(dòng)床內(nèi)密集顆粒處理成可壓縮流體,建立了單相可壓縮流μ(I)模型,對(duì)多個(gè)算例進(jìn)行了計(jì)算,并與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比.

1 基于改進(jìn)μ(I)模型的數(shù)學(xué)模型

1.1 控制方程

對(duì)于多相流系統(tǒng),采用多相流方法顆粒相的控制方程的張量表示如下

對(duì)于本文所研究的移動(dòng)床,沒(méi)有外部強(qiáng)制流入的氣流,床內(nèi)顆粒密集流動(dòng),顆粒密度遠(yuǎn)大于氣體密度,空隙中的空氣對(duì)顆粒的流動(dòng)影響很小,故氣相壓力p可忽略.顆粒固相擬壓力ps,一般認(rèn)為由碰撞動(dòng)力學(xué)壓力和摩擦壓力組成,確定固相擬壓力ps是顆粒擬流體模擬的關(guān)鍵點(diǎn)與難點(diǎn).移動(dòng)床內(nèi)顆粒密集,對(duì)于本文的研究對(duì)象,假設(shè)摩擦壓力pf起主導(dǎo)作用,只考慮摩擦壓力的影響,并采用Johnson 和Jackson[17](J &J)的公式描述摩擦壓力

對(duì)于移動(dòng)床內(nèi)密集顆粒流的計(jì)算,一般常取n=2,m=5,Fr=0.1αs[17],αs,min和αs,max分別代表顆粒最小和最大體積分?jǐn)?shù).

在多相流模型中,αsρs是顆粒的堆積密度,由于顆粒的體積分?jǐn)?shù)αs是變化的,故顆粒的堆積密度也在不斷變化.在對(duì)于本文研究的移動(dòng)床,忽略氣相對(duì)密度的影響,直接將αsρs視為擬流體密度

由式(3)和式(4),建立了顆粒流擬流體的密度ρ與顆粒體積分?jǐn)?shù)αs和顆粒擬壓力ps關(guān)系式,移動(dòng)床內(nèi)顆粒體積分?jǐn)?shù)變化,顆粒流密度及擬壓力隨之變化,顆粒流動(dòng)可表征為單相可壓縮流動(dòng).質(zhì)量力只考慮重力,控制方程為

1.2 本構(gòu)關(guān)系

應(yīng)用牛頓形式的張量表達(dá)描述移動(dòng)床內(nèi)摩擦應(yīng)力和應(yīng)變的關(guān)系

右側(cè)第一項(xiàng)為正應(yīng)力,在前文已有表述,應(yīng)用Johnson 和Jackson[17]的公式描述.

右側(cè)第二項(xiàng) τij切應(yīng)力.在μ(I)中,它被描述為摩擦系數(shù)和壓力的乘積,在本模型中,體系的壓力只考慮固相壓ps

μ為擬流體的摩擦系數(shù),由μ(I)關(guān)系確定.認(rèn)為μ為無(wú)量綱慣性數(shù)I值函數(shù)

慣性數(shù)表征是微觀時(shí)間尺度和宏觀時(shí)間尺度的比值.d/(ps/ρs)0.5為微觀時(shí)間尺度,它表示顆粒在壓力ps下落入大小為d中所需的時(shí)間,表征微觀重排的時(shí)間尺度.宏觀時(shí)間尺度1/γ表征宏觀平均變形時(shí)間尺度.I較小,一般認(rèn)為對(duì)應(yīng)準(zhǔn)靜態(tài)狀態(tài),I較大,進(jìn)入稠密慣性區(qū),顆粒之間主要作用為摩擦,I進(jìn)一步增大,則可能進(jìn)入快速流動(dòng),顆粒間相互作用以碰撞為主.需說(shuō)明的是,對(duì)于本文所研究的移動(dòng)床工況,床內(nèi)顆粒流動(dòng)低速而密集,假設(shè)顆粒壓力主要是摩擦壓力,不考慮顆粒的碰撞作用.

1.3 壁面邊界條件

黏性流體的無(wú)滑移的邊界條件不能用于描述顆粒流在壁面處的運(yùn)動(dòng),根據(jù)Coulomb 摩擦定律,壁面最大靜摩擦力為壓力和壁面摩擦系數(shù)的乘積.當(dāng)壁面剪切力小于壁面最大靜摩擦力時(shí),使用壁面剪切力;當(dāng)壁面剪切力大于壁面最大靜摩擦力時(shí),使用最大靜摩擦力.為便于編程計(jì)算,本文從Coulomb 摩擦定律出發(fā),類(lèi)比顆粒流內(nèi)部μ(I)模型,將壁面視作一個(gè)大顆粒,顆粒在壁面所受剪切力為

將壁面摩擦系數(shù)μw和慣性數(shù)I相關(guān)聯(lián)起來(lái),得到一個(gè)可小范圍變化的壁面摩擦系數(shù)

1.4 參數(shù)正則化及計(jì)算方法

在μ(I)模型中,由于擬流體壓力和黏度高度耦合的復(fù)雜關(guān)系,給數(shù)值求解帶來(lái)了困難,為保證數(shù)值求解的穩(wěn)定性和收斂性,近年來(lái)不少學(xué)者采用參數(shù)正則化方法[32-35]來(lái)解決此問(wèn)題.本文也對(duì)黏度和擬流體壓力采用正則化技術(shù).

將式(8)與式(9)之間引用黏度進(jìn)行過(guò)渡

將式(9)代入式(13)并整理并進(jìn)行正則化處理

式中右端第一項(xiàng)決定顆粒流的屈服應(yīng)力,第二項(xiàng)表示非牛頓黏性的貢獻(xiàn),λr為遠(yuǎn)小于1 的正則化參數(shù),本文計(jì)算中取λr=1 × 10-4.

對(duì)模型方程求解時(shí),必須保證黏度為正,所以擬流體壓力ps需要始終為正值,故同樣對(duì)壓力進(jìn)行正則化處理

上述模型與計(jì)算方法基于Fluent 軟件進(jìn)行,通過(guò)用戶(hù)自定義函數(shù)編程改寫(xiě)?zhàn)ざ群兔芏纫约斑吔鐥l件對(duì)模型進(jìn)行實(shí)現(xiàn).應(yīng)用壓力基求解器,采用壓力速度的耦合算法.

2 3 種典型散料在縮口料倉(cāng)內(nèi)的流動(dòng)

采用上述改進(jìn)的顆粒流單相可壓縮流μ(I)模型及控制方程,以表面光滑的球形玻璃珠、表面粗糙的球形剛玉球以及表面粗糙而形狀不規(guī)則的粗沙為研究對(duì)象,模擬研究了3 種典型散料在薄層移動(dòng)床縮口料倉(cāng)內(nèi)的出流流動(dòng)過(guò)程.

移動(dòng)床縮口料倉(cāng)結(jié)構(gòu)如圖1 所示,為便于實(shí)驗(yàn)測(cè)量,料倉(cāng)前后方向的厚度為30 mm,整個(gè)料倉(cāng)高500 mm,寬240 mm,上部為豎直段,顆粒在重力作用下向下流動(dòng);縮口段角度為60°,縮口斜面長(zhǎng)度為230 mm,出口寬度為10 mm.本課題組對(duì)不同散料在該裝置內(nèi)的流動(dòng)特性進(jìn)行了實(shí)驗(yàn)研究[36],采用PIV 測(cè)量了玻璃珠、剛玉球和粗沙在料倉(cāng)中速度.本文針對(duì)實(shí)驗(yàn)工況進(jìn)行了模型,初始取床內(nèi)顆粒體積分?jǐn)?shù)為0.5,出口根據(jù)實(shí)驗(yàn)測(cè)量采用速度出口,上部采用壓力入口條件,對(duì)于本文的計(jì)算條件,取入口固相壓力為400 Pa.計(jì)算中網(wǎng)格的劃分采用非結(jié)構(gòu)網(wǎng)格,網(wǎng)格數(shù)約18 萬(wàn),收斂條件為相對(duì)誤差小于0.5%.與采用多相流算法相比,由于避免了氣相流場(chǎng)和兩相耦合計(jì)算,計(jì)算時(shí)間大幅度減小.

圖1 料倉(cāng)模型Fig.1 Sketch of the silo

需說(shuō)明的是,料倉(cāng)前后方向的厚度如上文所述為30 mm,一般實(shí)驗(yàn)研究認(rèn)為這樣的薄板結(jié)構(gòu)可以忽略厚度方向的變化,稱(chēng)此種料倉(cāng)為近二維料倉(cāng),模擬計(jì)算中則不能避免前后壁面的影響,為三維計(jì)算,厚度方向上的網(wǎng)格數(shù)為14 個(gè),計(jì)算結(jié)果近前后壁面流速略低,其他位置變化不大.考慮到擬流體模擬中壁面效應(yīng)的誤差,本文在后續(xù)分析和與實(shí)驗(yàn)對(duì)比中所取數(shù)據(jù)為取距前壁第2 層網(wǎng)格計(jì)算結(jié)果.

μ(I)模型中計(jì)算參數(shù)與實(shí)際散料的實(shí)際物性密切相關(guān),本文所采用的參數(shù)具體如表1 所示,其中密度ρs、粒徑d參考實(shí)驗(yàn)測(cè)量值取值,I0和I0,w根據(jù)文獻(xiàn)[21]取為0.279.另需說(shuō)明的是,對(duì)于本文研究的薄層移動(dòng)床縮口料倉(cāng),實(shí)驗(yàn)測(cè)量床內(nèi)堆積狀態(tài)下顆粒的體積分?jǐn)?shù)在0.5 左右,為預(yù)報(bào)床內(nèi)體積分?jǐn)?shù)的變化,本文在計(jì)算中取αs,min和αs,max分別為0.45 和0.55.對(duì)于 μ (I) 模型中待確定的參數(shù)μs,μ2,μw,s,μw,2,基于實(shí)測(cè)物性,初步確定各個(gè)參數(shù)的取值范圍,進(jìn)行多組對(duì)比計(jì)算,和PIV 實(shí)驗(yàn)結(jié)果[36]相對(duì)照,最后確定了在改進(jìn)的 μ (I) 擬流體模型中玻璃珠、剛玉球和粗沙的μs,μ2,μw,s,μw,2參數(shù),具體如表1 所示,玻璃珠的內(nèi)摩擦系數(shù)以及與壁面的摩擦系數(shù)最低,剛玉球次之,粗沙最大.3 種散料的 μ (I) 模型參數(shù)的取值規(guī)律與實(shí)測(cè)顆粒的堆積角、壁面堆積角[36]的變化規(guī)律一致.

表1 模擬參數(shù)Table 1 Simulation parameters

圖2 是計(jì)算所得3 種散料的速度云圖.3 種散料速度分布差異十分明顯,玻璃珠由于內(nèi)部?jī)?nèi)摩擦系數(shù)以及與壁面的摩擦系數(shù)很小,在料倉(cāng)內(nèi)整體流動(dòng)最均勻,只在拐角近壁處出現(xiàn)明顯的低速區(qū),在其他區(qū)域玻璃珠都有較大流速.由于出口位于料倉(cāng)中心,在水平方向上,玻璃珠中心流速大,兩側(cè)較低,壁面處顆粒仍有較高流速;在豎直方向上,由于縮口和重力加速作用,玻璃珠在中心處的流速不斷增大,出口處達(dá)到最大值.剛玉球的內(nèi)部?jī)?nèi)摩擦系數(shù)以及與壁面的摩擦系數(shù)都較玻璃珠大,拐角附近以及豎直段兩側(cè)壁面處出現(xiàn)低速區(qū),流場(chǎng)的不均勻性增大.粗沙流場(chǎng)的不均勻性較剛玉球進(jìn)一步增大,近側(cè)壁有大范圍的顆粒滯止區(qū),顆粒主要在中心流道內(nèi)流動(dòng),先進(jìn)入料倉(cāng)的兩側(cè)顆粒需在中心顆粒的高速剪切攜帶下進(jìn)入中心流道才有可能流出料倉(cāng).大量的研究表明,不同顆粒在縮口料倉(cāng)中的流型一般可分為整體流和漏斗流,顆粒流動(dòng)性越好就越趨向于整體流.上述計(jì)算所得的流場(chǎng)分布中,玻璃珠的流動(dòng)屬典型整體流,粗沙則呈現(xiàn)明顯的漏斗流特征.這一結(jié)果也和文獻(xiàn)[36]中3 種散料的實(shí)驗(yàn)結(jié)果相一致.圖2 的結(jié)果表明,本文所采用的改進(jìn)的μ(I)模型及計(jì)算方法能有效揭示不同物性散料在縮口料倉(cāng)內(nèi)的流動(dòng)特性,對(duì)整體流和漏斗流的流動(dòng)特性都能較好表征.

圖2 不同散料的的速度分布Fig.2 Velocity distribution of different granular materials

Dai 等[36]采用PIV 測(cè)量了玻璃珠、剛玉球和粗沙在圖1 所示A-A,B-B,C-C3 個(gè)截面上的速度分布,其中B-B截面為縮口開(kāi)始位置,A-A和C-C分布位于其上、下100 mm.圖3 所示為3 個(gè)水平面上豎直方向上的速度模擬值與實(shí)驗(yàn)值的對(duì)比,散點(diǎn)是文獻(xiàn)中的實(shí)驗(yàn)數(shù)據(jù),實(shí)線是模擬結(jié)果.各個(gè)截面上模擬結(jié)果與實(shí)驗(yàn)整體吻合較好,但局部仍有一定誤差,如對(duì)于剛玉球,上部豎直段近壁處模擬速度低于實(shí)驗(yàn)值,粗砂中心流道的寬度低于實(shí)驗(yàn)值,出口最大略高于實(shí)驗(yàn)值.造成誤差的原因一方面是由于實(shí)驗(yàn)本身有一定誤差,另一方面,顆粒擬流體的假設(shè)也不可避免會(huì)犧牲局部流動(dòng)細(xì)節(jié)的預(yù)報(bào).

圖3 速度實(shí)驗(yàn)值與模擬的對(duì)比Fig.3 Velocity comparison between the experiment and simulation

圖4 所示為3 種散料的速度應(yīng)變率云圖,從速度應(yīng)變率的角度也明顯體現(xiàn)3 種散料內(nèi)部相互作用的不同.玻璃珠在料倉(cāng)的上部豎直段趨于整體流,顆粒之間的摩擦剪切作用較弱,在料倉(cāng)的上部豎直段部分基本沒(méi)有高應(yīng)變率剪切帶存在,下部縮口段受到拐角作用,形成低速區(qū),產(chǎn)生局部的應(yīng)變剪切帶.剛玉球和粗砂由于顆粒內(nèi)摩擦系數(shù)以及顆粒與壁面的摩擦系數(shù)大,近壁處顆粒滯流,中心主流區(qū)流速高,流場(chǎng)內(nèi)產(chǎn)生明顯的速度高應(yīng)變率剪切帶,剛玉球的剪切帶更窄且更靠近壁面.高應(yīng)變率的剪切帶將剛玉球和粗砂在料倉(cāng)上部豎直段內(nèi)的流動(dòng)分為顆粒低速區(qū)和顆粒主流區(qū).在底部縮口出流段,3 種散料都出現(xiàn)明顯的高應(yīng)變率區(qū)域,這是由于料倉(cāng)出口截面變小、顆粒流速變化大所引起的.

圖4 不同散料的速度應(yīng)變率分布Fig.4 Strain rate distribution of different granular materials

圖5 所示為3 種散料的慣性數(shù)分布云圖.由式(10)可知局部慣性數(shù)主要取決于局部速度應(yīng)變率,因此不同散料的慣性數(shù)分布特點(diǎn)與速度應(yīng)變率分布基本類(lèi)似,剪切應(yīng)變高的地方慣性數(shù)大.玻璃珠在料倉(cāng)上部豎直段慣性數(shù)較小,都在0.001 以下,剛玉球和粗砂則呈現(xiàn)兩條高慣性數(shù)帶狀區(qū)域,在高慣性數(shù)帶狀區(qū)域以及料倉(cāng)出口附近顆粒流的慣性數(shù)I≥ 0.001,其他區(qū)域顆粒流的慣性數(shù)都小于0.001.

圖5 不同散料的慣性數(shù)分布 (0～ 0.001)Fig.5 Inertial number distribution of different granular materials (0～ 0.001)

玻璃珠、剛玉球和粗沙在縮口料倉(cāng)內(nèi)流動(dòng)的慣性數(shù)范圍分別是10-7～ 0.23,10-7～ 0.5 和10-7～ 0.33.由于變化范圍大且大部分都在0.001 以下,圖5 的慣性數(shù)云圖的顯示范圍為0～ 0.001,圖6 所示為玻璃珠慣性數(shù)顯示范圍為0～ 0.1 的分布云圖,剛玉球和粗沙的云圖與玻璃珠基本相同.對(duì)于本文所計(jì)算的薄層移動(dòng)床縮口料倉(cāng)內(nèi)的出流流動(dòng),僅出口很小范圍內(nèi)顆粒流的慣性數(shù)大于0.1.有關(guān)顆粒流分區(qū)的定量劃分,現(xiàn)有文獻(xiàn)尚無(wú)統(tǒng)一的標(biāo)準(zhǔn)[13-15],對(duì)于恢復(fù)系數(shù)較大的顆粒,文獻(xiàn)[14] 認(rèn)為慣性數(shù)小于0.001時(shí)為準(zhǔn)靜態(tài)流動(dòng),在0.001 到0.1 之間是稠密慣性區(qū),大于0.1 是快速區(qū),若按此劃分,圖5 中紅色區(qū)域大部分為稠密慣性流動(dòng)區(qū)域,圖6 中出口的紅色小部分區(qū)域?yàn)榭焖賲^(qū),對(duì)于此區(qū)域本文的模型忽略顆粒碰撞造成的計(jì)算誤差可能較大,但如圖6 所示,此部分占比極小,本模型對(duì)絕大部分區(qū)域的模擬是適用的.

圖6 玻璃珠慣性數(shù)分布(0～ 0.1)Fig.6 Inertial number distribution of glass beads (0～ 0.1)

改進(jìn)的μ(I)擬流體模型的優(yōu)點(diǎn)在于可以對(duì)移動(dòng)床內(nèi)顆粒的體積分?jǐn)?shù)進(jìn)行預(yù)報(bào).圖7 所示為計(jì)算所得3 種物料的體積分?jǐn)?shù)在料倉(cāng)中的分布.玻璃的體積分?jǐn)?shù)在床內(nèi)變化較小,變化范圍分別是0.495～0.508,上部豎直段水平方向上變化較小,下部流動(dòng)通道變窄,顆粒體積分?jǐn)?shù)變化增大,拐角低速區(qū)顆粒流速低、體積分?jǐn)?shù)最大,到出口段,由于顆粒在重力作用下快速流出,顆粒體積分?jǐn)?shù)達(dá)到最小值.玻璃珠在下部縮口段水平方向上變化相對(duì)平緩,豎直方向上逐漸變小,出口處迅速降為最小值,這一趨勢(shì)是合理的,與文獻(xiàn)[37]定性相符.剛玉球和粗沙由于存在近壁低速區(qū)和流速較大主流區(qū),當(dāng)流速較大時(shí),顆粒的堆積也必然較為松散,所以主流區(qū)的體積分?jǐn)?shù)較低,滯流低速區(qū)體積分?jǐn)?shù)則相對(duì)較大,計(jì)算域內(nèi)剛玉球局部顆粒體積分?jǐn)?shù)最大達(dá)0.510,粗沙更可高達(dá)0.515,到出口處剛玉球和粗沙局部顆粒體積分?jǐn)?shù)最小達(dá)0.490.這可以從力鏈的角度來(lái)理解,在料倉(cāng)的滯流低速區(qū),顆粒之間相互作用主要通過(guò)力鏈傳播,最后會(huì)匯聚到壁面,壁面會(huì)承受很大一部分顆粒的重力,同時(shí)力鏈會(huì)使得壁面附近的顆粒群更為“壓實(shí)”,所以壁面附近的體積分?jǐn)?shù)更大.另外,物料在出口處體積分?jǐn)?shù)都為最小值,這一趨勢(shì)也是合理的.

圖7 不同散料的體積分?jǐn)?shù)分布Fig.7 Volume fraction distribution of different granular materials

3 移動(dòng)床顆粒繞管流動(dòng)特性研究

顆粒繞管移動(dòng)床在工業(yè)領(lǐng)域應(yīng)用廣泛,管道的布置對(duì)床內(nèi)顆粒流動(dòng)造成很大的影響,本文首先對(duì)移動(dòng)床單管繞流特性進(jìn)行了模擬研究,旨在分析改進(jìn)的μ(I)擬流體模型對(duì)存在單管影響的情況下能否有效描述.如圖8 所示,管道直徑30 mm,管道中心位于料倉(cāng)中心軸線上,離上部入口距離為205 mm.所研究的顆粒為玻璃珠,選用的模擬參數(shù)如表1 所示.由于圓管的阻擋,出口流速為低于上節(jié)出流工況,計(jì)算中取為20 mm/s,入口壓力邊界條件的設(shè)定方法與上節(jié)相同.

圖8 顆粒繞流流動(dòng)模型Fig.8 Sketch of granular flow around a pipe

移動(dòng)床單管繞流的速度和顆粒擬壓力云圖如圖9 和圖10 所示,圓管位于管道中心,對(duì)中心流動(dòng)顆粒的流動(dòng)起阻礙作用,圓管上部豎直段速度分布均勻,到管道上方速度變小,壓力增大,管道正下方則速度小,壓力小,管道兩側(cè)顆粒流速大.與一般黏性流體不同,除管正上方和正下方很小的區(qū)域,近壁區(qū)域的顆粒流速不為零,兩側(cè)顆粒以較大的速度流過(guò)圓管壁面,管道下方也沒(méi)有出現(xiàn)一般黏性流體的負(fù)壓回流區(qū)域,上述流動(dòng)特性與文獻(xiàn)[26]采用單相常密度μ(I)模型所得結(jié)果一致,也與文獻(xiàn)[36]實(shí)驗(yàn)觀測(cè)相符,說(shuō)明采用本文的模擬與方法能有效預(yù)測(cè)顆粒流移動(dòng)床單管繞流流動(dòng)的特征.

圖9 顆粒流繞管速度分布Fig.9 Velocity distribution of granular flow around a pipe

圖10 顆粒流繞管固相擬壓力分布Fig.10 Particle pressure distribution of granular flow around a pipe

圖11 所示移動(dòng)床單管繞流流場(chǎng)的體積分?jǐn)?shù)分布,單管上方由于顆粒流速低,體積分?jǐn)?shù)較大,局部最大體積分?jǐn)?shù)為0.510,下方體積分?jǐn)?shù)較小,管下方低體積分?jǐn)?shù)的區(qū)域與低固相壓力的區(qū)域相一致,正下方局部最小顆粒體積分?jǐn)?shù)為0.463.顆粒繞流經(jīng)過(guò)單管后進(jìn)入下部縮口段,顆粒流速增加,體積分?jǐn)?shù)減小,出口處局部顆粒體積分?jǐn)?shù)最小為0.461.

圖11 顆粒流繞管顆粒體積分?jǐn)?shù)分布Fig.11 Particle volume fraction distribution of granular flow around a pipe

顆粒流繞管體積分?jǐn)?shù)的變化近年來(lái)有不少研究報(bào)道,部分實(shí)驗(yàn)還觀測(cè)到了管下方的空隙區(qū),管下方空隙區(qū)的存在依顆粒材料和出口出流控制方式不同而變化,但目前普遍認(rèn)同顆粒在管下方體積分?jǐn)?shù)減小,本文的模型對(duì)管下方低體積分?jǐn)?shù)這一現(xiàn)象能較好模擬.但是,擬流體模擬由于忽略單個(gè)顆粒的大小,對(duì)于由于壁面效應(yīng)引起的貼壁處顆粒體積分?jǐn)?shù)變小無(wú)法精確預(yù)報(bào).

圖12 所示為顆粒流繞管速度應(yīng)變率云圖分布,管道的附近存在復(fù)雜變化的速度變化.管上部低速區(qū)外緣的顆粒受到主流區(qū)顆粒高速流動(dòng)的影響,出現(xiàn)高速度應(yīng)變率區(qū)域,圓管壁面對(duì)兩側(cè)顆粒的影響小,兩側(cè)顆粒保持高速流動(dòng),兩側(cè)對(duì)稱(chēng)出現(xiàn)兩個(gè)低速度應(yīng)變率區(qū)域,這部分高速流動(dòng)的顆粒又對(duì)管下部外緣區(qū)的顆粒造成剪切,因此管下方兩側(cè)形成很小的兩個(gè)高速度應(yīng)變率區(qū)域,但管道正下方的顆粒無(wú)法受到相鄰顆粒的裹挾,相應(yīng)地出現(xiàn)一個(gè)小的速度應(yīng)變率低值區(qū)域.

圖12 顆粒流繞管速度應(yīng)變率分布Fig.12 Strain rate distribution of granular flow around a pipe

顆粒流繞流經(jīng)過(guò)管道后逐步回到低速度應(yīng)變率水平,相應(yīng)的速度云圖趨于均勻.另外圖12 也表明了單管對(duì)周?chē)鲃?dòng)的影響情況,單管的存在使繞管兩側(cè)顆粒流速增大,對(duì)下部縮口段拐角處低速顆粒的攜帶剪切作用增大,使得拐角處的兩道剪切帶的位置更向壁面接近.

圖13 所示為移動(dòng)床單管繞流流場(chǎng)的慣性數(shù)分布,結(jié)合圖13(a)顯示范圍為0～ 0.001 和圖13 (b)顯示范圍為0～ 0.1 的慣性數(shù)分布云圖,按文獻(xiàn)[14]的劃分,上部豎直段大部分區(qū)域顆粒處于準(zhǔn)靜態(tài);繞管附近局部區(qū)域和下部縮口段速度應(yīng)變率增大,顆粒摩擦作用增大,大部分處于稠密慣性區(qū);出口處有很小區(qū)域內(nèi)慣性數(shù)大于0.1,處于顆?？焖倭鲄^(qū).雖然與圖6 相比,慣性數(shù)大于0.1 的區(qū)域變大,但在全流場(chǎng)占比仍然很小,本模型對(duì)絕大部分區(qū)域的模擬是適用的.

圖13 顆粒流繞管慣性數(shù)分布Fig.13 Inertia number distribution of granular flow around a pipe

4 結(jié)論

本文基于顆粒流單相μ(I)模型,采用J &J 摩擦壓力公式,補(bǔ)充局部顆粒體積分?jǐn)?shù)與顆粒局部壓力和局部顆粒流密度的關(guān)系式,將移動(dòng)床內(nèi)密集顆粒處理成可壓縮流體,建立了顆粒流單相可壓縮流μ(I)模型,并建立了顆粒流-壁面摩擦條件,在計(jì)算中對(duì)顆粒流黏度和壓力進(jìn)行正則化處理.

采用上述模型與算法對(duì)3 種典型散料在移動(dòng)床縮口料倉(cāng)內(nèi)的流動(dòng)進(jìn)行模擬,與實(shí)驗(yàn)對(duì)比,得到了玻璃珠、剛玉球和粗沙3 種散料的μ(I)模型參數(shù),并對(duì)3 種散料在料倉(cāng)內(nèi)的流動(dòng)特性進(jìn)行了分析,得到了3 種散料在料倉(cāng)內(nèi)的速度、應(yīng)變率和體積分?jǐn)?shù)等詳細(xì)分布,對(duì)于本文的計(jì)算工況,玻璃珠、剛玉球和粗沙在縮口料倉(cāng)內(nèi)流動(dòng)的顆粒體積分?jǐn)?shù)的變化范圍分別是0.495～ 0.508,0.490～ 0.510 和0.490～ 0.515,絕大部分區(qū)域流動(dòng)慣性數(shù)小于0.1,計(jì)算結(jié)果對(duì)玻璃珠在料倉(cāng)內(nèi)的整體流和粗沙的漏斗流特性能較好預(yù)報(bào),模擬結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果定性符合較好.

進(jìn)而對(duì)移動(dòng)床顆粒單管繞流流動(dòng)進(jìn)行了模擬,模擬結(jié)果詳細(xì)揭示了繞管周?chē)念w粒速度、應(yīng)變率和體積分?jǐn)?shù)分布.單管上方顆粒流速低、體積分?jǐn)?shù)較大,下方體積分?jǐn)?shù)小,計(jì)算結(jié)果合理,與文獻(xiàn)中揭示的顆粒繞管流動(dòng)特性相符.

上述計(jì)算結(jié)果表明,對(duì)于本文所研究的沒(méi)有外部強(qiáng)制氣流、顆粒在薄層移動(dòng)床內(nèi)低速密集流動(dòng)情況,顆粒體積分?jǐn)?shù)變化最大范圍為0.461～ 0.510,絕大部分區(qū)域流動(dòng)慣性數(shù)小于0.1,本文所提出的單相可壓縮流μ(I)模型與算法是可行的,計(jì)算量較顆粒離散元方法和多相流方法大幅度減小.當(dāng)然,基于擬流體假設(shè)的方法對(duì)由于壁面效應(yīng)所引起的局部顆粒體積分?jǐn)?shù)變化的預(yù)測(cè)、固相壓力模型、算法的收斂性等多方面還有待進(jìn)一步發(fā)展.