噴動床最小噴動速度的計算顆粒流體力學數值模擬

金海安， 王東祥,2， 俞建峰,2

(1.江南大學 機械工程學院， 江蘇 無錫 214122；2.江南大學 江蘇省食品先進制造裝備技術重點實驗室， 江蘇 無錫 214122)

噴動床具有結構簡單，床內氣固接觸和顆粒混合特性良好等優點，已廣泛應用于食品、化工、能源和輕工等領域[1-5]。噴動床有著特定的操作條件，當噴口表觀氣速達到最小表觀噴動速度(以下均稱最小噴動速度)，氣體才可以突破床層形成噴動。最小噴動速度對噴動床的工藝設計與優化具有重要意義。Mathur等[6]提出的粗顆粒最小噴動速度經驗關聯式應用最為廣泛。Monazam等[7]測量了普通平底噴動床的最小噴動速度，并基于Mathur-Gishler方程，建立了新的最小噴動速度關聯式。任立波等[8]耦合CFD-DEM算法，研究了錐形噴動床內的氣固流動行為，獲得了穩定流態化階段顆粒的時均速度分布規律。Duarte等[9]采用了雙歐拉模型對錐型噴動床內的流體動力學特性進行了研究，其最小噴動速度模擬結果與文獻試驗值吻合較好。近些年，通過CPFD方法對顆粒相進行雙重處理，節省了大量計算成本和時間，特別適合于復雜的多相流系統氣固流動行為的模擬研究，在工業級流化裝置模擬中逐漸成為主流[10-12]。課題組基于CPFD方法，建立二維錐形噴動床數值模型，數值模擬不同結構參數和操作參數床內的氣固流動特性；分析噴口寬度、噴動氣速以及顆粒特性對氣固流動結構和最小噴動速度的影響，建立二維錐形噴動床最小噴動速度預測模型。

1 CPFD方法

CPFD結合了歐拉和拉格朗日方法的優點，將顆粒相既視為連續介質，也視為離散單元，按歐拉方法以流體網格上的應力梯度代替顆粒應力梯度，再插值到離散顆粒體上，并對顆粒相的其他屬性使用拉格朗日法進行計算。CPFD將多個具有相同的物質屬性、物理運動及化學變化的真實顆粒打包成一個計算粒子[13]，計算速度快，并且可以保證全局和局部物質輸運守恒，可以更為快速準確地模擬復雜多相流系統的氣固流動行為。

1.1 控制方程

對于氣相，采用Navier-Stokes方程計算，湍流模型選擇大渦模型(LES)，流體連續性方程和動量方程如下[14]：

(1)

(2)

式中：uf為流體速度,m/s；θf為流體體積分數；ρf為流體密度,kg/m3；p為流體壓強,Pa；F為流體與顆粒相間單位體積的動量交換率；g為重力加速度,m/s2。

對于顆粒相，假設每個顆粒的質量在時間上是恒定的(顆粒之間或流體之間沒有質量傳遞)，但顆粒可具有一系列尺寸和密度，通過求解顆粒分布函數Φ(xp,up,ρp,Ωp,t)的輸運方程來描述顆粒相的運動。其中xp為粒子位置；up為粒子速度，m/s；ρp為粒子密度，kg/m3；Ωp為粒子體積，m3。輸送方程和顆粒的加速度[15]為：

(3)

(4)

式中：A是計算粒子加速度，m/s2；θp是顆粒體積分數；τp是顆粒的法向應力，Pa；Dp是曳力系數。

1.2 顆粒碰撞模型

CPFD中將多個真實顆粒處理成一個計算顆粒，對于真實顆粒間的碰撞無法逐一進行計算，通過引入碰撞模型來表征真實顆粒間的碰撞，將流體網格中的應力梯度插值到相應位置的顆粒上進行計算。顆粒碰撞時法向應力τp[16]為

(5)

式中：Ps為常數，Pa；θcp為顆粒堆積時的體積分數；β為常數，推薦值2～5；ε取值為10-7。

1.3 曳力模型

多相流中氣固間的相互作用力通過曳力模型表征，常用的有Wen-Yu模型[17]、Ergun模型[18]和Gidaspow模型[19]等，其中Wen-Yu模型適用于稀相顆粒流，Ergun模型更多的適用于密相區域，Gidaspow模型結合了Wen-Yu和Ergun 2個模型的特點，其適用的濃度范圍更為寬泛。噴動床內氣固流動呈現稀相和密相混合狀態，課題組選擇Gidaspow模型計算曳力系數Dp：

Dp=D1,θp<0.75θcp;

(6)

0.75θcp<θp<0.85θcp; (7)

Dp=D2, 0.85θp<θp。

(8)

式中：D1和D2分別為Wen-Yu模型和Ergun模型的曳力系數。

2 模型參數

2.1 數值模型

如圖1(a)所示，噴動床床高為800 mm，床體寬度為200 mm，厚度20 mm。噴動氣入口底部是夾角為60°的V形結構，入口厚度為20 mm。模擬中入口厚度不變，入口寬度分別調整為20，15和10 mm以改變噴動氣的進口尺寸。采用笛卡爾網格生成方法進行網格劃分，如圖1(b)所示。如果在網格中形成非常小的切割單元不能在模型中包含相關數量的計算粒子，則必須被移除或與附近較大的單元合并來維持模型穩定性。研究采用的顆粒直徑為1.3～3.2 mm，床體寬度是床體厚度的10倍，床體厚度大于粒徑的5倍，滿足二維噴動床條件[20]，其呈現出的流動特性能較真實反應實際的氣固流動。

采用網格數量為10 000，15 000和20 000進行網格無關性驗證以檢驗合理的網格數量范圍。從表1中可以看出，當網格數量超過15 000后，床體高度為50，100和200 mm處壓力監測點的數值不再出現明顯改變，即可保證計算結果的可靠性。為了更準確地模擬噴動床內顆粒的流動形態，底部錐形區域的網格進行了加密，有效網格數量為16 251。

表1 網格無關性驗證

2.2 模擬參數

以空氣作為入射氣體，所選固體物料為綠豆、小米和3種不同直徑的玻璃珠，顆粒參數如表2所示。表3所示為模擬時相關參數的設置情況。采用升速法結合壓降曲線得到最小噴動速度并對顆粒不同流動形態進行觀察。

表2 顆粒物性參數

表3 模擬參數

3 結果與討論

3.1 實驗驗證

根據Liu等[21]的實驗數據對初始床高分別為50，75和100 mm的二維噴動床最小噴動速度進行模擬計算，填料為直徑2 mm，密度2 380 kg/m3的玻璃珠，通過逐漸增加噴動氣速引起的床層壓降變化和具體床層流動形態來確定最小噴動速度。圖2為表觀氣速Us=1.58 m/s時穩定噴動狀態下二維噴動床的實驗和模擬對照圖，從圖中可以看出兩者在噴動結構上有著高度一致性。環隙區顆粒不斷流入噴動區，最終在噴動區出口形成噴泉。模擬結果中噴泉區中顆粒相對較少的主要原因是CPFD方法是將多個真實顆?！按虬背梢粋€計算顆粒進行計算。

3.2 噴動結構

圖4為Us=2.3Ums時噴動過程的模擬情況。t=0.0 s時，向噴動床內通入高速氣體來推動物料上升；t=0.2 s時，床層中部形成一段大氣泡狀空腔，尺寸接近床層截面積，類似騰涌射流流動結構。主要原因是氣速突然增大，推動物料整體上升，床層表面高度明顯升高，氣體需要運動距離增大，向四周滲透的氣體增加，造成空腔范圍增大。氣體繼續上升至突破床層后將攜帶的顆粒向四周拋灑，逐漸形成穩定的噴動區和噴泉區?？梢钥吹江h隙區顆粒緩慢向下運動至底部后被卷吸入噴動區隨著入口高速氣流上升，進行新的循環。

圖5是噴動氣速逐漸增大時床內流動結構的變化。當噴動氣速小于最小噴動氣速時，氣體無法穿透床層，顆粒只能在床層內部形成射流區內進行循環。區別于圖4中t=0.2 s時產生的大氣泡，圖5中Us=0.8Ums時形成的內部射流區較窄，在內部射流區頂端形成氣泡繼續上升并不斷聚攏至床層表面破裂。隨著噴動氣速繼續增加，氣流突破床層后并形成噴泉。從圖中可以明顯看出噴泉的高度隨著噴動氣速增加而增加，當氣體速度足夠大時，大量顆粒被吹起布滿整個床體，沒有明顯噴泉區界線，噴動床開始向氣力輸送床轉變。

3.3 最小噴動速度的影響因素

3.3.1初始床高的影響

圖6為噴口寬度為20 mm的條件下，初始床高和床寬的比值與最小噴動速度的關系。從圖中可以看出在噴口寬度和床寬不變的條件下，不同填充材料的最小噴動速度都隨著初始床高的增加而增加，并且在初始床高超過錐形體高度后，兩者呈現出良好的正比關系。產生這種現象的原因主要有2個：一是因為初始床高的增加即初始填充物料的增加，氣體在上升過程中需要吹起更多的顆粒，導致阻力增大，就需要更高的氣速來克服阻力產生噴動；二是初始床高的增加使氣體需要上升的距離變長，更多的氣體在突破床層的過程中由噴動區向四周環隙區逸散，導致需要更多氣體來產生噴動，使最小噴動速度增大。

3.3.2顆粒直徑的影響

圖7為顆粒直徑與最小噴動速度的關系，在噴口寬度為20 mm的條件下，選用了1.3，1.8和2.3 mm的3種直徑的玻璃珠作為填充物料。從圖中可以看出最小噴動速度隨著顆粒直徑的增大而迅速增大。主要原因是由于同種密度的材料，顆粒直徑增大，顆粒的表面積和體積比減小，導致相同質量的顆粒噴動所需氣體的速度增加[22]。隨著顆粒直徑的增大，顆粒的終端速度相應增大也是最小噴動速度增大的一個原因。結合圖6初始床高與最小噴動速度的關系可以看出同樣的比例系數下顆粒粒徑對最小噴動速度的影響遠大于初始床高的影響。

3.3.3噴口寬度的影響

圖8為噴口寬度與最小噴動速度的關系。在H0/Dt為1.4的條件下，噴口寬度為10，15和20 mm時，3種填充物料的最小噴動速度變化呈現出相同趨勢。從圖中可以看出隨著噴口寬度增大，最小噴動速度逐漸增大。一方面是隨著噴口寬度的增加，噴動區的平均寬度也逐漸增大，氣體在上升過程中會卷吸更多的顆粒，需要更高的氣速來克服其產生的阻力。另一方面，由于噴口寬度的增加，氣體向環隙區的擴散量也在增加，導致氣體攜帶的動能損失增大，需要更多的氣體來彌補動能損失，導致最小噴動速度增大。

通過上述分析得到最小噴動速度受到初始床高、顆粒直徑、顆粒密度和噴口寬度等影響。根據唐鳳翔等[23],[24]54對最小噴動速度關聯式的研究，可以寫成下面的參數群形式：

(9)

結合不同操作條件下的模擬結果進行多元線性回歸分析，確定系數a～e，得到最小噴動速度關聯式為

(10)

圖9(a)將最小噴動速度關聯式的計算值與模擬值進行比較，從圖中可以看出計算值與模擬值具有很高的吻合性；相關系數為0.98，標準誤差為9.8%，吻合較好，相對誤差在12.5%以內，該關聯式對此類噴動床最小噴動速度的預測具有較高的準確性。圖9(b)中將得到的最小噴動速度的關聯式與其他學者的關聯式進行對比，從圖中可以發現在相同的條件下，所有關聯式的預測值呈現相同趨勢；相對誤差大部分控制在30%以內，與Uemaki和Wang等[25-26]得到關聯式預測值非常吻合，說明研究得到的關聯式對于最小噴動速度的預測具有一定的可靠性。Zhong[24]52的關聯式中考慮了流化氣的影響，而Choi等[27]填充顆粒的粒徑比差異較大，噴動床床體并非二維結構，導致了與關聯式的預測值誤差較大。

4 結論

課題組采用CPFD方法建立了二維噴動床數理模型，對最小噴動速度和氣固流動結構進行數值模擬，考慮了操作參數和結構參數的改變對最小噴動速度和流動結構轉變的影響，并通過對比文獻實驗數據驗證了模型準確性，得到如下結論：

1) 隨著噴動氣速增加，床內將逐步呈現內部射流、鼓泡射流、噴動、穩定噴動和氣力輸送等流動結構，穩定噴動流動結構的發生需要噴動氣速維持在一定范圍內。

2) 最小噴動速度隨著初始床高、顆粒粒徑和噴口寬度增大而增大，通過增大噴動速度可克服顆粒產生的阻力和彌補氣體攜帶動能損失；相對于噴口寬度，靜止床高和顆粒直徑對最小噴動速度影響較為顯著。

3) 研究得到的二維錐形噴動床最小噴動速度關聯式與已有二維經驗關聯式吻合較好，能夠準確預測二維噴動床最小噴動速度，對于三維噴動床則存在一定誤差。