基于雙流體模型的噴動-流化床不同流動形態的模擬研究☆

趙俊楠,王會寧,戈朝強,劉國棟

(1.哈爾濱工業大學 能源科學與工程學院，黑龍江 哈爾濱 150001;2.北京航天動力研究所,北京 100076)

0 引言

在噴動床(Spouted Bed)中，高速噴動的氣體通過噴嘴噴入床體，形成高速射流，射流攜帶顆粒噴動至一定高度，顆粒受到重力作用，在達到最高處之后向床體兩側自由下落，進入環狀回流區，受到高速射流的夾帶再次向上噴動形成循環。由于大尺寸顆粒可以在噴動床中取得較好的流化效果，因此被廣泛應用于農作物干燥、生物制藥等重要領域。但是，這種床型在床底兩側會形成不流化不流動的“死區”。噴動-流化床(Spout-fluid Bed)的出現，顯著提升了床體兩側的流化程度，原先“死區”內部的顆粒被充分流化起來，得到了更為廣泛的應用。目前，噴動-流化床技術在煤氣化過程中扮演著重要角色，成為煤氣化技術的主要爐型之一。除此以外，在農業生產、生物制藥等多個領域，噴動-流化床技術也有相當廣泛的應用。

近年來，國內外的專家學者已經對噴動-流化床進行了深入而廣泛的研究。但對其流動結構以及動力學參數的掌握依然不夠系統，在工業放大和實際應用等過程中仍然有一定困難。因此，本文的研究目的在于運用雙流體模型研究噴動-流化床內不同流動結構下的氣固兩相流動特點，探討兩相的濃度、速度以及脈動速度的分布規律。

長期以來，很多學者采用不同方法對氣固兩相噴動流化系統進行了研究[1-3]。沈來宏等[2]采用實驗方法研究了噴動氣速和流化氣速對噴動-流化床內流化特性的影響。Kuipers等人[3]采用雙流體的方法模擬了三維流化床提升管，得到了氣固兩相的流動特點。李乾軍等[4]采用雙流體模型研究了噴動-流化床的相關流化特性，模擬結果中，床層上部出現清晰的噴泉狀結構，床層上下有明顯的分層現象，床層兩側有明顯的環核結構。唐楠等[5]采用雙流體模型研究了不同流化速度對噴動-流化床流化特性的影響。肖睿等[6]采用雙流體模型研究了加壓導向式噴動流化床，研究發現顆粒由于受到壓力影響在噴動區加速上升，環形區域的氣體速度在卷吸顆粒時變大，而空隙率則有下降趨勢。Mathur[7]以及唐鳳翔等[8]研究了噴動床內的最小噴動速度，其研究表明，最小噴動速度與床層高度成正比，與床體直徑成反比。

考慮到研究對象的特性，以及氣固兩相之間的相互作用力，本文采用歐拉雙流體模型，通過顆粒動力學和摩擦應力模型建立稠密氣-固兩相模擬的本構方程，以研究不同噴動速度、流化速度對噴動-流化床內流動結構、流速以及顆粒脈動的影響。

1 數值計算模型

1.1 控制方程

質量方程，對于氣相

(1)

對于固相

(2)

動量守恒方程，對于氣相

(3)

對于固相

(4)

1.2 相間相互作用

本文采用Gidspow[9]曳力模型，它結合了Wen &Yu[10]和Ergun[11]模型各自的優點，其表達式如下所示

(5)

(6)

其中，單顆粒曳力系數Cd取決于顆粒Re數

(7)

(8)

1.3 應力模型

為了封閉固相動量方程(方程(4))，需要給出固相壓力ps和應力τs的表達式。根據Johnson & Jackson[12]的研究，固相壓力ps和應力τs都由運動-碰撞和摩擦兩部分構成

ps=pkc+pf

(9)

τs=τkc+τf

(10)

運動-碰撞部分的固相應力張量為

(11)

式中μb——體積粘度;

μkc——剪切粘度。

摩擦部分的固相應力張量為

τf=2μfDs

(12)

式中μf——摩擦粘度;

Ds——應變張量的偏量。

(13)

氣相動量方程(方程(3))中的氣相應力張量為

(14)

其中，氣相應變量為

(15)

1.3.1 顆粒動理學

運動-碰撞部分的固相壓力pkc與固相粘度μkc通常由顆粒動理學[13]進行封閉

pkc=[1+2εs(1+e)g0]ρsεsθ

(16)

(17)

(18)

(19)

其中，γs為碰撞耗散

(20)

κs為湍動能擴散率

(21)

本文采用Bagnold[14]的徑向分布函數

(22)

1.3.2 摩擦應力模型

在固相高濃度區域，顆粒之間的相互作用不再是瞬時碰撞，而是長時間的摩擦。因此，采用顆粒動理學來封閉固相應力明顯不合適，需要引入摩擦應力模型[15]來描述顆粒之間的長時間擠壓與摩擦

(23)

(24)

式中φ——內摩擦角。

2 計算條件

2.1 模型結構及邊界條件

圖 1為本文所模擬的噴動-流化床結構示意圖，單位為mm。噴口尺寸為22 mm×12 mm，噴動氣流由此噴入流化床內部，底部周圍區域為流化風入口，通過調節噴動風速度和流化風速度可以獲得不同的流動工況，圖1所示的結構尺寸單位均為mm。頂部為壓力出口，與大氣壓力相等。其余部分設為壁面，其中氣體為無滑移壁面邊界條件，顆粒為部分滑移壁面邊界條件。

圖1 模型結構示意圖

網格尺寸為顆粒尺寸的2倍。初始床層高度為0.228 m，其他模擬參數詳見表1。

表1 模擬參數

本文的模擬結果要與Link等人[1]的實驗結果進行對比，因此工況的選取與Link等人[1]一致，包含了5種流動結構，具體參數見表2。

表2 工況參數

注：ubg為流化風速度，usp噴動風速度。

3 計算結果及分析

3.1 流動結構對比分析

(1)工況A：ubg=1.8 m/s；usp=85 m/s

此條件下的流型為Spouting with aeration，數值模擬顆粒濃度云圖和Link 等[1]實驗快照對比如圖2所示。從圖中可以看出，在床內中部存在明顯的氣體通道，氣體通道會穿透整個床層。由于噴動風速度遠遠大于兩側流化風速度，大部分顆粒堆積在床內兩側，只有氣體通道中及其附近的顆粒才被流化，隨著噴口氣體的噴出而呈噴泉狀運動，屬于典型的噴動流化運動模式。

圖2 工況A顆粒濃度云圖與實驗快照對比

(2)工況B：ubg=2.0 m/s；usp=76 m/s

此條件下的流型為Intermediate，與工況一相比，噴動風速度與流化風速度的比值由47.22降至38。數值模擬顆粒濃度云圖和Link等[1]實驗快照對比見圖3。從圖中可以看出，雖然降低了噴動風速度，但在床內中部仍然存在明顯的氣體通道，空隙率較大。由顆粒速度矢量圖可以看出，床內顆粒速度較大的位置在顆粒氣體通道內及其周圍，且顆粒主要向上運動，而在兩側處顆粒有較小的向下運動趨勢，在底部流向氣體噴道，形成對稱性較好的循環運動。

圖3 工況B顆粒濃度云圖與實驗快照對比

(3)工況C：ubg=2.2 m/s；usp=66 m/s

此條件下的流型為Spout-fluidization，噴動風速度與流化風速度的比值降至30，此時床層上部的顆粒被流化并開始緩慢移動。數值模擬顆粒濃度云圖和Link等[1]實驗快照對比見圖4。由圖4可以看出，在床層表面下方氣體噴道中可以觀察到兩側顆粒形成頸部結構。可能是由于大量顆粒流化噴射后，自由下落至環狀回流區對噴道產生擠壓作用，造成顆粒周期性地阻塞氣體通道。

圖4 工況C顆粒濃度云圖與實驗快照對比

(4)工況D：ubg=2.8 m/s；usp=37 m/s

此條件下的流型為Jet in fluidized bed；Link等實驗快照見5(a)，數值模擬顆粒濃度云圖見圖5(b)。從圖中可以看出，這時幾乎所有顆粒都在運動，而且在噴口兩側的環狀區域中連續形成氣泡。氣泡在床層底部產生，隨著高度增加逐漸長大，上升至一定高度，氣泡破裂顆粒先被揚起，然后在重力的作用下自由下落，氣固混合運動較為劇烈。但此時，氣泡尺寸較小，其形成和傳播的時間較短。這種流動結構類似于鼓泡床中的流化特性。但此時仍然還存在一個氣體通道，該通道周期性地被大量的顆粒阻塞，且通常發生在床層的低位。

圖5(a) 工況D,Link等實驗快照圖[1]

圖5(b) 工況D數值模擬顆粒濃度云圖

(5)工況E：ubg=3.5 m/s；usp=3.5 m/s

此條件下的流型為Slugging bed，此時所有顆粒都在運動并形成了不均勻結構-顆粒聚團。圖6展示了顆粒聚團形成的不同階段， Link 等的實驗快照見圖6(a)，本文數值模擬顆粒濃度云圖見圖6(b)。從圖中可以看出，顆粒形成聚團并周期性的向上運動，然后在重力作用下下落，數值模擬結果與實驗結果吻合。除了顆粒聚團之外，在床層底部還形成了較小的氣泡，沿著床層高度逐漸長大，上升至自由表面后發生破裂。

圖6(a) 工況E,Link等實驗快照圖[1]

圖6(b) 工況E數值模擬顆粒濃度云圖

3.2 流動參數分析

3.2.1 顆粒速度分布

圖7表示在0.05 m、0.1 m、0.15 m和0.2 m高度處，不同工況下顆粒的徑向速度分量。在0.15 m以上的區域，工況A至工況D的顆粒由噴口向兩側噴出，下落并回流至環狀區，這一區域顆粒不再繼續大幅度上升，而是向兩側擴散，屬于噴泉區。在0.1 m以下區域，從噴泉區落下的顆粒進入環狀回流區，由于噴口處的顆粒被噴動風帶走，環狀區域的顆粒會流向噴口形成循環。環狀回流區的顆粒徑向速度剛好與噴泉區相反。工況E的徑向速度均為正值，說明顆粒是形成聚團整體移動。

圖7 不同高度處的顆粒徑向速度分量

圖8表示在0.05 m、0.1 m、0.15 m和0.2 m高度處，不同工況的顆粒軸向速度分量。如圖所示，工況A至工況D的軸向速度曲線在中軸處存在一個波峰。速度的最大值與觀察高度成反比，這是因為觀察高度越高，噴動風的速度衰減得越多，且顆粒自身重力作用也會使顆粒速度減小。由噴口往外顆粒軸向速度逐漸減小，在靠近壁面處為負值，這是因為壁面附近的顆粒在回落顆粒的擠壓下而向下運動。曲線關于中軸對稱分布兩個波谷，這是因為四周環空區下落的顆粒在壁面處由于壁面的摩擦作用，速度會有所減小，下降速度的最大值大約在離壁面0.02 m的位置。

圖8 不同高度處的顆粒軸向速度分量

工況E條件下，軸向速度曲線較平緩沒有十分陡峭的波峰波谷，這是因為顆粒聚團整體運動，在y軸方向有相近的速度。

3.2.2 顆粒擬溫度分布

顆粒擬溫度用來表征顆粒脈動能，數值越大表明顆粒脈動越劇烈。圖9表示在0.05 m、0.1 m、0.15 m和0.2 m高度處，不同工況下的顆粒擬溫度的分布曲線。如圖所示，在0.1 m以下區域，工況A至工況D顆粒擬溫度曲線均在中軸處出現一個陡峭的波峰，遠離噴口中心的顆粒擬溫度較低，趨近于0。這是因為中心處為噴口，高速射流的噴動風夾帶著顆粒劇烈運動，因此脈動劇烈。四周的流化風速度較低，不會使顆粒發生劇烈脈動。此區域工況E的顆粒擬溫度曲線十分平緩且接近于0，說明該條件下，這一區域的顆粒脈動都十分微弱。

圖9 不同高度處顆粒擬溫度分布

在0.15 m以上的區域，工況A至工況D顆粒擬溫度曲線在中軸處仍存在一個波峰。與0.1 m以下區域相比，波峰的最大值由原來的0.3減小到0.17 m2/s2左右；且波峰的徑向跨度變寬；靠近壁面處的顆粒擬溫度升高，使得曲線出現了兩個波谷。這是因為越往上，噴動風速度衰減越多，顆粒的脈動強度也隨之降低；噴動區上升的顆粒至一定高度后回落，與下方上升的顆粒、氣流以及壁面發生碰撞，使得顆粒的湍動加劇。工況E的顆粒擬溫度沿徑向仍較平緩，曲線的中間段仍趨近于0，只是在靠近壁面處略有上升，有可能是氣流沿壁面上升穿過顆粒聚團導致的。

4 結論

本文采用Euler-Euler雙流體模型，建立了布風板為平板型的噴動流化床計算模型，對噴動流化床內氣固兩相流動進行了數值模擬計算，調節流化風速度和噴動風速度模擬計算了5個工況，分別對應于噴動流化床的5種流型，得到了床內顆粒速度、濃度以及顆粒擬溫度等參數分布。

模擬得出的顆粒濃度分布與J.M.Link等[1]的實驗快照對比效果理想，清晰直觀的反映出噴動流化床濃度分布規律；尤其對于Spout-fluidization流型在噴道中部會形成頸部結構，以及Slugging bed流型會形成不均勻結構——顆粒聚團，模擬計算濃度云圖能夠清晰觀察到這些結構。

在噴射區，顆粒主要沿豎直方向向上運動，速度的最大值與噴動風速度成正比；在噴泉區，顆粒主要沿水平方向向壁面擴散，而豎直方向的速度已大幅降低；在環隙區，顆粒在重力及噴動風影響下沿豎直方向下運動，但下降速度并不大，且在壁面附近由于壁面摩擦作用，會出現速度降低的“壁面效應”。

在噴動流化床的低位，噴射區的顆粒擬溫度最高，顆粒擬溫度的最大值也與噴動風速度成正比，在環隙區顆粒擬溫度幾乎等于0；在噴動流化床的高位，噴泉區中心的顆粒溫度較高，但是遠小于噴動區中心的顆粒擬溫度，與噴動流化床低位的環隙區不同，噴泉區的壁面附近顆粒擬溫度略有回升，這是向四周擴散的回落顆粒速度脈動引起的。

符號表

拉丁字母

t——時間

pg,ps——氣相，固相壓力

Cd——單顆粒曳力系數

dp——顆粒直徑

Re——顆粒雷諾數

pkc——運動-碰撞顆粒壓力

Ds——固相應變張量的偏量

Dg——氣相應變張量

I——單位張量

e——恢復系數

g0——徑向分布函數

ks——固相湍動能擴散率

pf——固相摩擦壓力

希臘字母

ρg,ρs——氣相，固相密度

εg,εs——氣相，固相濃度

εs,max——填充極限

τg,τs——氣相，固相應力張量

τkc——固相運動-碰撞應力

τf——固相摩擦應力

μb——固相體積粘度

μkc——固相運動-碰撞粘度

μf——固相摩擦粘度

β——曳力系數

φ——內摩擦角

θ——顆粒擬溫度

γs——湍動能碰撞耗散率