李慧君， 王業(yè)庫， 張久意

（華北電力大學能源動力與機械工程學院，河北保定 071003）

鼓泡流化床由于具有氣固兩相接觸與混合充分、溫度分布均勻和傳熱傳質(zhì)效率高等優(yōu)點，已被廣泛應用于工業(yè)生產(chǎn)中［1-4］。流化床內(nèi)的氣固流動本質(zhì)上是動態(tài)的，通過實驗研究流化床費用較高，條件苛刻，很難獲取有效的實驗數(shù)據(jù)。隨著計算機模擬技術的快速發(fā)展，計算流體力學（CFD）被廣泛用于研究鼓泡流化床的氣固流動特性。

在氣固流動模型中，歐拉-歐拉雙流體模型（TFM）因計算量較小且準確性較高而被廣泛使用［5］。氣固兩相間的作用力非常復雜，一般情況下，相較于曳力，其他作用力可忽略不計，因此僅考慮曳力作用，進而曳力系數(shù)的選擇對于流化床的模擬精度影響很大［6］。TFM包含常用的曳力模型，即Wen Yu模型、Ergun模型和Gidaspow模型等均勻模型，其中前兩者均有其合適的顆粒體積分數(shù)范圍；而Gidaspow模型雖然沒有具體顆粒體積分數(shù)范圍限制，但沒有考慮氣泡和顆粒團簇等非均勻結(jié)構。氣泡和顆粒團簇的特性很大程度上會影響流化床的流體動力學和運行效率［7-10］，因此引出介尺度結(jié)構的亞網(wǎng)格曳力模型［11］。李靜海等［12］在考慮顆粒團簇的基礎上提出了能量最小多尺度（EMMS）模型，該模型以及以此為基礎的其他模型均獲得了很高的精度［13-15］。Shi等［16］基于鼓泡流化床提出EMMSBubbling模型，發(fā)現(xiàn)該模型能精準預測出床內(nèi)的顆粒分布狀態(tài)以及顆粒軸向和徑向速度的變化。

EMMS模型雖然模擬精度較高，但受到顆粒本身物性參數(shù)的影響，且沒有普適性，很難應用于實際工程。筆者根據(jù)不同曳力模型的特點，在不同顆粒體積分數(shù)區(qū)間選擇不同的曳力模型。為克服新型曳力模型不連續(xù)性的問題，在顆粒體積分數(shù)分界點引入光滑函數(shù)使其連續(xù)，在此基礎上耦合TFM，對使用典型B組顆粒的二維鼓泡流化床進行模擬，并將其模擬結(jié)果與Gidaspow模型的模擬結(jié)果以及實驗值進行對比，選取最適合B組顆粒的修正因子。

1 數(shù)學模型

顆粒整體的運動表現(xiàn)出類似流體的運動規(guī)律，可以把顆粒相視為具有相同顆粒性質(zhì)的一種擬流體［17］，采用TFM［18］結(jié)合顆粒動力學（KTGF）［19］對鼓泡流化床進行模擬，通過求解氣固兩相的各項連續(xù)性方程、動量守恒方程以及用來封閉方程的本構方程，獲得模擬結(jié)果。

1.1 控制方程

氣相質(zhì)量方程［17］為：

式中：φ為體積分數(shù)；ρ為密度，kg／m3；u為平均速度，m／s；mgs為氣固間的質(zhì)量源項，kg／（m3·s）；t為時間，s；下標g、s分別表示氣相和固相（即顆粒）。

氣相動量方程為：

式中：τ為剪切應力，Pa；p為壓力，Pa；β為曳力，kg／（m3·s）；g為重力加速度，m／s2。

氣相能量方程為：

式中：λ為導熱系數(shù)，W／（m·K）；H為比焓，J／kg；T為溫度，K；hgs為氣相與固相間的傳熱系數(shù)，W／（m2·K）。

固相質(zhì)量方程為：

固相動量方程為：

固相能量方程為：

式中：α為空隙率。

式中：I為單位張量；kΘs為能量的擴散系數(shù)；γΘs為擬顆粒溫度因非彈性碰撞引起的能量耗散系數(shù)；φgs為顆粒隨機波動動能。

1.2 本構方程

各相剪切應力τi

［18］為：

式中：V為瞬時速度，m／s；ξ為體積黏度，Pa／s，對于氣相而言，ξ為0 Pa／s；μ為剪切黏度，Pa／s；下標i可表示g或s。

固相壓力ps為：

式中：go為徑向分布函數(shù)；ess為顆粒彈性碰撞后的恢復系數(shù)。

固相體積黏度ξs為：

式中：d為直徑，m。

徑向分布函數(shù)為：

式中：Φ為內(nèi)摩擦角，（°）；I2D為偏應力張量第二部變量，N／m2。

能量的擴散系數(shù)為：

擬顆粒溫度因非彈性碰撞引起的能量耗散系數(shù)為：

顆粒隨機波動動能為：

考慮到顆粒的聚團作用，將鼓泡流化床分為密相（φs＞0.4）、過渡相（0.1≤φs≤0.4）和稀相（φs＜0.1）。密相采用Ergun模型，稀相采用Wen Yu模型［19］，過渡相采用經(jīng)過介尺度修正的Mckeen模型［11］。3種曳力模型計算式分別為：

式中：βErgun、βWen Yu和βMckeen分別為Ergun模型、Wen Yu模型和Mckeen模型的曳力，kg／（m3·s）；Re為雷諾數(shù)；C為修正因子，對于大顆粒常取0.5～0.99。

上述3種曳力模型在固相體積分數(shù)分界點處并不連續(xù)，因此引入光滑函數(shù)ψ1和ψ［19］2，使其成為連續(xù)的曳力模型。

利用式（23）和式（24）連接3種曳力模型，重新構建新型曳力βgs的數(shù)學模型，其計算式為：

鼓泡流化床常用B組顆粒作為流化顆粒，因此新型曳力的修飾因子分別選取0.7、0.8和0.9。同組顆粒具有相似的物理性質(zhì)。通過將TFM與新型曳力模型和Gidaspow曳力模型耦合，對比不同曳力模型模擬結(jié)果與實驗值的誤差，得到B組顆粒的修正因子。

2 幾何模型

根據(jù)文獻［20］中的實驗進行參數(shù)設定，以驗證新型曳力模型的正確性。該實驗將3個差壓傳感器置于距床層底部0.03 m、0.3 m和0.6 m的位置上，橫跨氣體分布器。在表觀氣速為0～0.8 m／s的條件下對總壓降和床層膨脹率進行監(jiān)測，當流化床流動達到穩(wěn)態(tài)后，記錄總壓降和床層膨脹率的測量值。使用反射式光纖探頭PC-4粉末空隙測量儀對床層空隙率進行測量，該探測器包含一束光纖，由于顆粒直徑比光纖束直徑小得多，光線被測量床層中的顆粒反射，從輸出電壓可獲得瞬時固相體積分數(shù)。該實驗在高為1.0 m、寬為0.28 m、厚度為0.025 m的二維有機玻璃柱中進行。球形玻璃微珠直徑為250～300μm，密度為2 500 kg／m3，屬于典型的B組顆粒，在環(huán)境條件下采用空氣進行流化。靜態(tài)床層高度為0.4 m，固相體積分數(shù)為0.6。采用二維簡化鼓泡流化床進行模擬，模擬時間為8 s。二維鼓泡流化床模型示意圖見圖1。入口邊界條件為速度入口，出口邊界條件為自由出口，時間步長為0.001 s，最大迭代次數(shù)為20，收斂準則為0.001。具體給定的模擬參數(shù)見表1。

圖1 二維床模型的示意圖及網(wǎng)格Fig.1 Schematic diagram of the two-dimensional bubbling fluidized bed model with its grid

表1 主要模擬參數(shù)和條件Tab.1 Main simulation parameters and conditions

3 結(jié)果分析

3.1 網(wǎng)格敏感性分析

選取不同網(wǎng)格數(shù)進行模擬，其中靜態(tài)床層區(qū)域的網(wǎng)格劃分與其他區(qū)域相同，并采用Gidaspow曳力模型進行評估。不同網(wǎng)格數(shù)下軸向顆粒體積分數(shù)分布見圖2。由圖2可知，3種網(wǎng)格數(shù)下軸向顆粒體積分數(shù)分布趨勢一致。網(wǎng)格數(shù)為2 800時，床層底部的軸向顆粒體積分數(shù)相比其他2種網(wǎng)格數(shù)較小，隨著網(wǎng)格數(shù)的增加，床層底部軸向顆粒體積分數(shù)逐漸增大。與網(wǎng)格數(shù)為2 800時相比，網(wǎng)格數(shù)為11 200時軸向顆粒體積分數(shù)的相對偏差約為23.86%；網(wǎng)格數(shù)為28 000時軸向顆粒體積分數(shù)的相對偏差約為29.17%?？紤]到模擬精度和計算時間，最終選用網(wǎng)格數(shù)為11 200。

圖2 不同網(wǎng)格數(shù)下軸向顆粒體積分數(shù)分布Fig.2 Distribution of axial particle volume fraction under different grid numbers

3.2 鼓泡流化床穩(wěn)定時域分析

在流態(tài)化初期，鼓泡流化床內(nèi)流動狀態(tài)不穩(wěn)定，不利于流態(tài)化分析，因此需排除初始流態(tài)的影響。圖3為采用Gidaspow模型和新型曳力模型獲得的瞬時壓差Δp1曲線。其中，瞬時壓差Δp1為距床層底部0.03 m與距床層底部0.3 m壓力監(jiān)測線的壓力差。由圖3可知，在3 s后不同模型的瞬時壓差基本穩(wěn)定波動，表明流化床內(nèi)顆粒已達到穩(wěn)定流態(tài)化。其中，壓差波動是由氣泡穿過床層時發(fā)生分裂、碰撞和合并造成的。

圖3 采用新型曳力模型和Gidaspow模型得到的瞬時壓差曲線Fig.3 Instantaneous pressure difference obtained by new drag model and Gidaspow model

3.3 床層時均壓差

采用Gidaspow模型和帶有不同修正因子的新型曳力模型模擬床層時均壓差Δp，并將其與實驗值進行比較，結(jié)果見圖4。其中，床層時均壓差Δp為距床層底部0.03 m與0.6 m處壓力監(jiān)測線的壓差。當表觀氣速較?。║g＜1.5Umf）時，模擬值與實驗值吻合度較差。其原因是在模擬過程中表觀氣速過小，顆粒沒有被充分流化，氣固兩相間作用力以摩擦力為主，導致沿氣體運行方向氣體壓差過大。在實驗過程中，顆粒并非均勻分布，如果有顆粒堵塞在進氣噴口附近，則進氣口壓力降低，導致床層時均壓差實驗值較小。當表觀氣速較大時，床層時均壓差模擬值與實驗值的吻合程度很高。表觀氣速為0.38 m／s時，修正因子為0.9的新型曳力模型的床層時均壓差模擬值與實驗值最接近，其誤差約為1.8%；當表觀速度達到0.46 m／s時，修正因子為0.8的新型曳力模型床層時均壓差的精度較Gidaspow模型提高了約3.1%；新型曳力模型模擬得到的床層時均壓差始終優(yōu)于Gidaspow模型，其中修正因子為0.8和0.9的新型曳力模型預測所得床層時均壓差精度相似，且優(yōu)于修正因子為0.7的新型曳力模型。

圖4 不同表觀氣速下床層時均壓差模擬值與實驗值的對比Fig.4 Comparison between simulated and experimental values of average pressure difference at different apparent gas velocities

3.4 床層膨脹率

采用床層膨脹率E［19］作為衡量新型曳力模型準確性的依據(jù)。

式中：H′為床層膨脹高度，m；H′0為靜止床層高度，m。

不同表觀氣速下床層膨脹率模擬值與實驗值的對比見圖5。由圖5可知，隨著表觀氣速的增大，不同曳力模型下床層膨脹率模擬值均逐漸提高。表觀氣速為0.1 m／s時，Gidaspow模型和帶有不同修正因子的新型曳力模型的床層膨脹率模擬值相差較小；表觀氣速為0.38 m／s時，修正因子為0.9的新型曳力模型的床層膨脹率模擬值相比Gidaspow模型提高了4.7%，在其他2種修正因子下新型曳力模型的床層膨脹率模擬值略大于Gidaspow模型；當表觀氣速為0.51 m／s時，修正因子為0.9的新型曳力模型的床層膨脹率模擬值相比Gidaspow模型提高了7.1%，修正因子為0.7和0.8時則分別提高1.8%和1.1%。隨著表觀氣速的不斷增大，新型曳力模型的床層膨脹率模擬精度不斷提升，其中修正因子為0.9的新型曳力模型預測所得床層膨脹率精度最高，修正因子為0.7和0.8的新型曳力模型的預測值接近。當φs為0.1～0.4時顆粒聚團嚴重，經(jīng)過介尺度修正的Mckeen模型相比Gidaspow模型模擬精度更高。

圖5 不同表觀氣速下床層膨脹率模擬值與實驗值的對比Fig.5 Comparison between simulated and experimental values of bed expansion rate at different apparent gas velocities

3.5 顆粒體積分數(shù)

圖6給出了Ug=0.38 m／s時瞬時顆粒體積分數(shù)云圖。由圖6可知，不同曳力模型模擬所得瞬時顆粒體積分數(shù)與實驗具有相似性。從圖6（a）可以看出，實驗中床層底部附近有小氣泡，當氣泡上升到頂部表面發(fā)生合并時，氣泡變大。氣泡的生長是由壁面效應以及與其他氣泡的相互作用導致的。曳力模型模擬與實驗觀測到的氣泡大小基本一致，但仍存在誤差。這是因為在實驗中氣體分布器的設計對其附近的射流穿透和流體流動特性有一定影響，而模擬過程中設置為均勻進氣，未考慮氣體分布器對氣泡的影響。

圖6 瞬時顆粒體積分數(shù)云圖Fig.6 Instantaneous particle volume fraction nephogram

從圖6還可以看出，不同曳力模型下靜態(tài)床層高度模擬值低于實驗值。其原因如下：（1）實驗時顆粒為非球形，導致實際曳力偏大；（2）TFM耗散相對更強，模擬的氣固結(jié)構不易維持；（3）實驗時顆粒直徑并非定值，而模擬時取為定值。

圖7給出了床層高度為0.2 m截面處不同表觀氣速下的時均顆粒體積分數(shù)變化。與Ug=0.38 m／s相比，當Ug=0.46 m／s時不同曳力模型的時均顆粒體積分數(shù)分布更加對稱，但其與實驗值之間的誤差較大。當Ug=0.38 m／s時，孔隙度曲線出現(xiàn)輕微的不對稱。這是由于床層中形成了一定的流型以及顆粒在床層上均勻分布的時間太短。當Ug=0.46 m／s時，在靠近床層中心處模擬和實驗得到的時均顆粒體積分數(shù)曲線均較為平坦，說明表觀氣速越大，流動越均勻。由圖7可知，新型曳力模型與Gidaspow模型得到的時均顆粒體積分數(shù)分布相似，其中修正因子為0.8的新型曳力模型的時均顆粒體積分數(shù)模擬值與實驗值誤差率最小，Gidaspow模型的誤差率最大。當Ug=0.38 m／s時，Gidaspow模型的最大誤差率約為21%；修正因子分別為0.8、0.7和0.9的新型曳力模型的最大誤差率分別約為11.1%、16%和17.8%；當Ug=0.46 m／s時，Gidaspow模型的最大誤差率約為33%，修正因子分別為0.8、0.7和0.9的新型曳力模型的最大誤差率分別約為18.1%、19.1%和19.3%。綜上，當Ug=0.38 m／s時，修正因子為0.8的新型曳力模型精度最高，比Gidaspow模型提高了約10%；當Ug=0.46 m／s時，修正因子為0.8的新型曳力模型精度最高，比Gidaspow模型提高了約15%。綜合床層時均壓差、床層膨脹率和顆粒體積分數(shù)，修正因子為0.8的新型曳力模型最適用于B組顆粒。

圖7 不同表觀氣速下的時均顆粒體積分數(shù)變化Fig.7 Variation of time average particle volume fraction at different apparent gas velocities

4 結(jié) 論

（1）隨著表觀氣速的增大，新型曳力模型得到的床層時均壓差模擬值逐漸接近實驗值。當表觀氣速增大至0.46 m／s時，修正因子為0.8的新型曳力模型精度比Gidaspow模型提高了3.1%。

（2）隨著表觀氣速的增大，新型曳力模型得到的床層膨脹率模擬精度優(yōu)于Gidaspow模型。表觀氣速為0.51 m／s時，修正因子為0.9的新型曳力模型的床層膨脹率模擬值與實驗值之間的誤差最小，相較于Gidaspow模型，其精度提高了7.1%。

（3）表觀氣速為0.38 m／s時，Gidaspow模型的時均顆粒體積分數(shù)最大誤差率約為21%，而修正因子為0.8的新型曳力模型的最大誤差率約為11.1%；當表觀氣速為0.46 m／s時，修正因子為0.8的新型曳力模型的時均顆粒體積分數(shù)分布曲線最優(yōu)，其精度較Gidaspow模型提高了約15%。

（4）新型曳力模型的預測精度均高于Gidaspow模型。修正因子為0.8和0.9的新型曳力模型精度接近，且優(yōu)于修正因子為0.7的新型曳力模型。修正因子為0.8的新型曳力模型最適用于B組顆粒。