基于介尺度曳力模型雙組分顆粒混合的CFD模擬

宋曉皎， 王 帥

(1. 山西能源學院能源與動力工程系， 山西晉中 030600； 2. 哈爾濱工業大學能源科學與工程學院， 黑龍江哈爾濱 150001)

流化床反應器被廣泛應用于催化裂解、 CO2捕集和煤粉顆粒燃燒等工業化生產中[1-3]。顆粒直徑和密度等物理特性的不同，會導致多組分氣固流化床中顆粒的混合分離現象的產生[4]，對氣固反應特性及燃料的轉化率產生直接影響，因此，對氣固流化床中雙組分顆粒分離特性的深入研究，將有利于獲取系統最佳操作條件，進而提高燃料利用率。

近年來，數值模擬成為研究多組分顆粒流化床系統的重要方法[5-7]。Feng等[8]利用離散顆粒模型對流化床內雙組分顆粒的分離和混合過程進行數值模擬， 并指出氣固相間的作用力引發了流態化和顆粒間分離現象的產生。 Owoyemi等[9]利用數值模擬的方法， 對三維流化床內金紅石顆粒的混合和分離過程進行研究， 并通過實驗證實模擬結果的正確性， 指出多組分顆粒混合系統中固-固相間作用力和固相剪切力需進行進一步研究。 Mazzie等[10]利用矩陣直接求解法， 分析了流化床內多組分顆粒混合現象， 并對數值擴散對模型預測的影響進行了討論。 Sharma等[11]討論了顆粒密度和直徑等物理特性對生物質和生物炭顆粒混合特性的影響， 結果顯示生物質顆粒密度對混合物的分離和混合特性影響最為顯著。

顆粒分離現象主要受流化床內氣泡移動和顆粒湍動的影響。基于麥克斯韋速度耗散假設，非均一顆粒擬溫度顆粒動力學理論被提出，并成功應用于雙組分顆粒的鼓泡流化床的數值研究中[12-13]。Iddir等[14]基于能量非均一性考慮，進一步發展、 完善了混合物顆粒動力學理論，并指出該理論能較好地應用于剪切流和反應流系統中。Chao等[15]考慮到顆粒間的摩擦影響，得到雙組分顆粒動力學理論，模擬分析顯示雙組分上升管中主要呈現出連續性流體的脈動現象。

雙組分氣固曳力模型的構建對于多組分氣固流化床系統的數值模擬具有十分重要的作用。Beetstra等[16]基于格子玻爾茲曼模擬的方法獲得應用于多組分系統的校正因子。Olaofe等[17]使用此校正因子對雙組分流化床內顆粒的混合和分離特性進行數值研究，發現源于直接數值模擬(direct numerical simulation, DNS)方法的曳力關系式可以獲得更好的模擬預測結果。Zhang等[18]應用離散顆粒方法，通過比較顆粒分離程度和氣泡生成頻率對多組分曳力模型的影響進行了研究；但此方法的應用受到顆粒直徑和系統尺寸等多方面的限制，因此，非常有必要發展一種考慮到床層非均一性、適用于大尺寸系統粗網格模擬的曳力模型[19-20]。

Shi等[21]提出了以氣泡為介尺度的能量最小介尺度曳力模型，并將其應用于鼓泡流化床的數值研究中。Lunge和Wang等[22-23]將曳力模型拓展應用于雙組分氣固流化床中，并且對雙組分混合系統中的顆粒分離混合現象以及反應特性進行研究，但忽略了其他顆粒相的影響。Zhou等[24]針對2類顆粒聚團現象對曳力模型進行改進，并將改進后的模型應用于循環流化床提升管的數值研究中，預測結果與實驗結果較吻合。然而，基于能量最小多尺度 (energy minimization multi-scale，EMMS)曳力模型、同時考慮乳化相中雙組分顆粒的相互影響及氣泡影響并應用于鼓泡流化床中的雙組分氣固曳力模型的相關研究較為少見。

本文中提出基于乳化相中顆粒相互影響并考慮以氣泡為介尺度的曳力模型，結合顆粒動力學理論和歐拉(Euler)多流體模型，對具有雙組分顆粒的三維鼓泡流化床中顆粒的混合和分離特性進行模擬研究，并通過實驗數據驗證雙組分氣固曳力模型的合理性，并進一步分析研究外界操作條件和顆粒物性參數對雙組分顆粒的分離和混合過程的影響。

1 數學模型

假定2種固相分別為單一直徑和密度的球形顆粒處于等溫流場內，建立基于歐拉多流體模型和顆粒動力學理論、 以氣泡為介尺度、 考慮乳化相中兩相顆粒相互影響的雙組分氣固曳力模型。在驗證雙組分氣固曳力模型準確性的基礎上，可利用此曳力模型對雙組分顆粒分離和混合過程進行數值模擬。

1.1 控制方程

在Euler多流體模型[23]中，固體顆粒和氣體被作為擬流體和連續性介質進行處理。其質量守恒方程為

(1)

(2)

式中：φg為氣相體積分數；ρg為氣相密度， kg/m3；Ug為氣相表觀速度， m/s；xi為三維空間中某一方向的坐標，i=1，2，3；m、n為反應物質編號；Ng為氣相反應物種類總數；Rgn為氣相間轉化的單位體積質量流率， kg/(m3·s)；φm為固相體積分數；ρm為固相密度， kg/m3；Umi為三維空間中某一方向的固相表觀速度，i=1，2，3，m/s；Nm為固相反應物種類總數；Rmn為固相間轉化的單位體積質量流率，kg/(m3·s)。

其動量守恒方程為

(3)

(4)

式中：Ugi和Ugj為三維空間中2個方向的氣相表觀速度，i=1，2，3，j=1，2，3，m/s；xj為三維空間中某一方向的坐標，j=1，2，3，m；Pg為氣相壓強，Pa；τgij為二維應力張量，Pa;gi為三維空間中某一方向重力加速度，i=1，2，3，m/s2；β為曳力系數；Umi和Umj為三維空間中2個方向的固相表觀速度，i=1，2，3，j=1，2，3，m/s；Pm為固相壓強，Pa；為向量微分算子符號。

其顆粒擬溫度方程為

(5)

式中：Θm為固相顆粒擬溫度，m2/s2；∏m為由于固相間滑移而產生的動能，J/(m3·s)；γ為因相間滑移而導致的能量碰撞耗散，J/(m3·s)；km為固相顆粒運動碰撞的熱傳導，J。

1.2 氣泡介尺度曳力模型

本文中考慮氣泡介尺度及乳化相中顆粒相的相互作用對鼓泡流化床內相間曳力的影響，在Shi等[21]提出的單組分氣固曳力模型的基礎上，構建適用于雙組分混合物流化床系統的氣固曳力模型。基于氣泡內無顆粒相存在的假設，得到乳化相的平均物性[25]為

ρe=ρp1φep1+ρp2φep2+ρgφeg，

(6)

μe=μg[1+2.5(1-φep)+10.05(1-φep)2+0.002 73 exp(16.6-16.6φep)]，

(7)

Ue=(ρgUge+ρp1Uep1+ρp2Uep2)/ρe，

(8)

式中：ρe為乳化相的平均密度，kg/m3；ρp1和ρp2分別為乳化相中顆粒p1和p2的密度，kg/m3；φep1和φep2為乳化相中顆粒p1和p2的體積分數；μe和μg分別為乳化相和氣體的動力黏度，N·s/m2；φep為乳化相顆粒總體積分數；Ue為乳化相表觀氣體速度，m/s；Uge、Uep1和Uep2為乳化相中氣相、顆粒p1和p2的表觀氣體速度，m/s。

通過引入加速度項來表征非穩定特性，氣泡相和乳化相中顆粒的力平衡公式為

(9)

(10)

式中：db為氣泡直徑，m；Cdb為傳統曳力模型的氣泡相有效曳力系數；Ueb為乳化相中氣泡表觀氣體速度，m/s；ab為氣泡相的加速度，m/s2；Cde1和Cde2分別為傳統曳力模型乳化相中顆粒p1和顆粒p2的有效曳力系數；ae為乳化相的加速度，m/s2。

為了表征氣泡相和乳化相的附加質量力和慣性的差異[26]，其差值為

(11)

式中：ρb為氣泡相的平均密度，kg/m3；δb為氣泡相體積分數；σ為局部顆粒體積分數濃度變化量。

基于氣體和固體體積分數的定義以及氣固質量守恒關系，可得到關系式：

(12)

Ug=Uge(1-δb)+δbUb，

(13)

(14)

式中：Usp1和Usp2分別為顆粒p1和顆粒p2的表觀滑移速度，m/s。

為封閉模型求解，設定輸送懸浮顆粒所需能量最小作為顆粒流動穩定準則[21]，即

(15)

式中：fb為氣泡中氣體體積與整個床層氣體體積分數比值。基于以上方程的求解，可獲得相關局部結構參數，進而得到雙組分顆粒各固體相的曳力系數表達式

(16)

(17)

式中：Usp1和Usp2為2種顆粒表觀滑移速度，m/s。基于此，為更好地表征考慮氣泡介尺度后曳力模型的變化情況，引入非均質系數Hd， 氣泡介尺度曳力系數計算式變為

(18)

Hd=βb/βw，

(19)

式中：Cd為單個顆粒的曳力系數；βb和βw分別為氣泡介尺度曳力系數和傳統WEN-YU曳力系數。

1.3 混合指數

為了更加深刻地分析多組分顆粒分離和混合特性，引用Owoyemi等[27]提出的混合系數定義式

M=wt/wo，

(20)

式中：wt為床層頂部(25%)下沉顆粒的平均質量分數；wo為整個床層中下沉顆粒的平均質量分數。

1.4 幾何模型及計算參數

參考Sun等[28]模擬試驗中的幾何模型，建立的三維流化床反應器模型的結構示意圖如圖1所示。為了研究顆粒的混合和分離特性，以砂礫顆粒S1和稻殼顆粒S2作為床料，2種顆粒的平均密度分別為2 600、 950.6 kg/m3，平均直徑分別為0.44、 1.54 mm；氣相采用無滑移邊界條件，固相采用Sinclair和Jackson邊界條件[29]。

圖1 三維流化床反應器模型的結構示意圖Fig.1 Schematic diagram of three dimensional fluidized bed reactor model

三維模擬計算的初始條件以及物料的材料屬性如表1和表2所示。

表1 三維模擬的初始條件

表2 材料屬性

2 結果與分析

2.1 非均質系數

為更直觀地呈現入口氣體表觀速度對兩相顆粒非均質系數的影響，在Sun等[28]選取入口氣體表觀速度為0.79 m/s實驗條件的基礎上，添加入口氣體表觀速度為1.50、 2.00 m/s這2種工況進行分析對比。入口氣體表觀速度、床層空隙率和顆粒體積分數對非均質系數的影響如圖2所示。由如圖2可知，非均質系數隨床層空隙率近似呈指數形式增長，不同顆粒相體積分數增長速率有所不同，這說明雙組分顆粒的非均質系數不僅僅取決于床層空隙率，還會受到不同粒徑的顆粒體積分數的影響，所以，利用單組分顆粒介尺度曳力模型對雙組分顆粒系統進行數值模擬是不合理的；除此之外，從圖2還可明顯觀察到，隨著氣體速度的增大，固相非均質系數呈逐漸減小趨勢，說明較高氣速時介尺度結構對曳力大小的影響會更顯著。

a)砂礫b)稻殼圖2 入口氣體表觀速度、床層空隙率和顆粒體積分數對非均質系數的影響Fig.2 Effects of inlet gas apparent velocity, bed voidage and particle volume fraction on heterogeneity coefficient

不同顆粒體積分數下非均質系數隨床層空隙率的變化如圖3所示。由圖3可知，兩相顆粒的非均質系數都隨著顆粒直徑的增大逐漸增大；隨著床層空隙率的增大，砂礫顆粒的非均質系數增長幅度逐漸減小，而稻殼顆粒直徑非均質系數增長情況則相反；較小直徑顆粒的體積分數變化對較大直徑顆粒的非均質系數變化趨勢的影響較大，并且直徑較小顆粒的非均質系數相對較小，說明較小直徑或質量的顆粒對相間曳力的非均質性影響程度更大，相間曳力與非均質結構的關聯性更強。

a)砂礫b)稻殼圖3 不同顆粒體積分數下非均質系數隨床層空隙率的變化Fig.3 Variation of heterogeneity coefficient with bed voidage under different particle volume fraction

部分流化床反應器并不總是處于常溫工作狀態，溫度不同必然會影響到氣相物理特性的變化，間接影響床層非均質特性及相間曳力。不同床層溫度條件下非均質系數隨床層空隙率的變化如圖4所示。由圖4可知，隨著床層溫度的升高，固相非均質系數逐漸減小，稻殼顆粒的非均質系數減小程度較小，這主要是由于床層溫度的增大直接導致氣體黏度增大，密度減小，使得相間曳力發生變化，導致床層溫度變化對較輕顆粒的非均質系數影響較大，而對較重顆粒的非均質系數影響較小。

a)砂礫b)稻殼圖4 不同床層溫度條件下非均質系數隨床層空隙率的變化Fig.4 Variation of heterogeneous coefficient with bed voidage at different bed temperatures

2.2 鼓泡流化床雙組分顆粒混合模擬結果

參照文獻 [28]中在入口氣體表觀速度分別為0.58、 0.79 m/s這2種工況下進行的雙組分顆粒混合實驗，利用氣泡為介尺度的雙組分氣固曳力模型和傳統Gidaspow曳力模型進行數值模擬，稻殼質量分數沿床層分布情況的模擬結果與實驗結果對比如圖5所示。由圖5可知，在2種入口氣體表觀速度條件下，床層頂部的模擬結果均略微小于實驗結果，但是整體上模擬結果與實驗結果吻合度較高；在床層上部氣泡介尺度曳力模型模擬結果比傳統Gidaspow曳力模型模擬結果更接近實驗結果，這主要是由于床層上部氣泡尺寸更大，對相間曳力的影響更加明顯，因此，以氣泡介尺度的雙組分氣固曳力模型能夠獲得更準確的預測結果，后續研究均采用該模型。

考慮到2種顆粒動力學特性差異，為更直觀地呈現入口氣體表觀速度對兩相顆粒混合和分離特性的影響，選取入口氣體表觀速度分別為0.31、 0.79 m/s， 稻殼顆粒直徑分別設為0.15、 0.44、 1.54 mm和床層溫度分別為298、 398、 498 K為后續數值模擬的初始條件。

a)氣體入口表觀速度為0.58 m/sb)入口氣體表觀速度為0.79 m/s圖5 稻殼質量分數沿床層分布情況模擬結果與實驗結果的對比Fig.5 Comparison of simulation results and experimental results of rice husk mass fraction distribution along bed

在t=7 s(t為進氣時間，以下同)時刻， 入口氣體表觀速度對瞬時顆粒濃度分布云圖的影響如圖6所示。 由圖6可知， 入口氣體表觀速度為0.31 m/s時， 床層內部的氣泡數量較多并且體積較小； 當入口氣體表觀速度增加到0.79 m/s時， 稻殼顆粒在床層上部的顆粒濃度相對高于床層底部的濃度， 但是兩相顆粒分離程度不明顯， 說明隨著入口氣體表觀速度的增大， 氣泡數量減少， 體積增大， 固相顆粒混合程度逐漸加大， 床層膨脹高度增大。 這是由于， 入口氣體表觀速度的增大使得氣泡快速聚集成大直徑氣泡， 到達稀密兩相分界面處時發生破裂， 對顆粒攜帶作用增強， 進而使得床內固相混合程度加大， 床層膨脹高度增加。

a)入口氣體表觀速度0.31 m/sb)入口氣體表觀速度0.79 m/s圖6入口氣體表觀速度對瞬時顆粒濃度分布云圖的影響Fig.6 Influence of inlet gas apparent velocity on instantaneous particle concentration distribution neogram

當入口氣體表觀速度為0.31 m/s時，稻殼顆粒直徑分別設為0.15、 0.44、 1.54 mm，在t=7 s時刻，稻殼顆粒直徑對瞬時顆粒濃度分布云圖的影響如圖7所示。由圖7可知，隨著稻殼顆粒直徑的減小，顆粒最小流化速度逐漸較小，使得兩相顆粒分離程度逐漸加大，雙組分顆粒的分離和混合還取決于顆粒的物性。

在不同入口氣體表觀速度條件下， 稻殼顆粒直徑對顆粒混合指數隨的影響如圖8所示。 由圖8可知， 隨著稻殼顆粒直徑的減小， 混合指數逐漸減小， 分離程度逐漸加大； 當入口氣體表觀速度減小時， 混合指數也隨之減小， 因此， 顆粒物性和外界操作條件的合理匹配將是實現雙組分顆粒高效分離的關鍵之舉。

a)dp2=0.15 mmb)dp2=0.44 mmc)dp2=1.54 mm圖7 稻殼顆粒直徑對瞬時顆粒濃度分布云圖的影響Fig.7 Effect of rice husk particle diameter on instantaneous particle concentration distribution neogram

圖8 稻殼顆粒直徑對顆粒混合指數的影響Fig.8 Effect of grain diameter of rice husk on grain mixing index

流化床反應器工作環境有時處于較高的溫度。 在入口氣體表觀速度為0.31 m/s時， 設定床層溫度分別為298、 398、 498 K， 在t=7 s時刻的床層溫度對瞬時顆粒濃度分布云圖的影響如圖9所示。 由圖9可知， 隨著床層溫度的升高， 氣體黏度增大， 顆粒相容易黏結， 進而導致兩相顆粒分離程度逐漸減弱。

床層溫度對顆粒混合指數的影響如圖10所示。由圖10可知，隨著床層溫度的升高，混合指數逐漸增大，分離程度逐漸減弱。 這主要是由于隨著床層溫度的升高， 氣體黏度增大， 密度減小， 使得氣固曳力系數減小，最終導致非均質系數減小； 兩相顆粒的非均質系數下降幅度不同，使得差值逐漸減小； 同時，氣體黏度增大，也使得顆粒較易黏結，不易分離，因此，除了通過改變顆粒物性和外界操作條件以外，還可以適當降低床層溫度來提高顆粒分離程度。

a)298 Kb)398 Kc)498 K圖9 床層溫度對瞬時顆粒濃度分布云圖的影響Fig.9 Effect of bed temperature on instantaneous particle concentration distribution nephogram

圖10 床層溫度對顆粒混合指數的影響Fig.10 Effect of bed temperature on particle mixing index

3 結論

以氣泡為介尺度建立了雙組分氣固曳力模型，在三維鼓泡流化床反應器中模擬了雙組分顆粒的分離和混合過程，得到以下結論：

1)非均質系數隨床層空隙率近似呈指數增長；固相顆粒直徑的減小、氣體入口表觀速度的增大以及床層溫度的升高都將會加強兩相間曳力的非均質性；相比于傳統的Gidaspow氣固曳力模型，雙組分氣固曳力模型的預測精度更高。

2)隨著入口氣體表觀速度的增大，固相顆粒混合程度逐漸加大，床層膨脹高度增大。隨著稻殼顆粒直徑和入口氣體表觀速度的減小，兩相顆粒分離程度加大，混合指數隨之減小；顆粒物性和外界操作條件的合理匹配有利于雙組分顆粒的分離；同時，隨著床層溫度的升高，兩相顆粒分離程度逐漸減弱，混合指數逐漸增大。

在今后的研究中，將進一步優化雙組分氣固曳力模型，并將其拓展應用到反應性流化床系統中，模擬分析流化床內的氣固反應和傳質特性。