基于Wray-Agarwal 湍流模型的鼓泡流化床氣固流動的模擬研究與分析

喬紅旺 宋釗毅 向 陽

(北京化工大學 化學工程學院, 北京 100029)

引 言

隨著現代工業的快速發展,人類社會對化石能源的過度依賴造成了一系列的生態環境問題。 在碳中和背景下,亟待開發清潔能源。 太陽能具有清潔、存量巨大等特點,是一種非常理想的清潔能源。 太陽能的開發利用需要太陽能電池板,其主要材料是高純度的多晶硅。 現階段多晶硅生產主要采用第三代西門子法,但硅烷流化床法具有結構簡單、能耗低等優點,是多晶硅產業未來發展的重要方向。 然而,流化床反應器內氣固兩相流體力學特征與氣相沉積反應動力學都很復雜,現有研究還遠遠不夠充分,需要更加深入的探討[1-2]。

近年來,隨著計算機技術的快速發展,計算流體力學(CFD)得到了更為廣泛的應用,在獲得實驗中難以測量的詳盡參數、闡明過程機理的同時,又大幅度降低了科研成本,安全性更高。

氣固流化床內多相流動是一種非常復雜的湍流流動,目前應用廣泛的k-ε和k-ω兩方程模型被用于反應器數值模擬研究。 Mehdizad 等[3]基于Shear Stress Transfer(SST)k-ω模型建立了一種新的數值模型,使其能夠在廣泛的顆粒直徑和顆粒密度范圍內研究氣固床中的最小流化速度。 Du等[4]選擇Standardk-ε湍流模型研究了多晶硅流化床反應器中多晶硅顆粒的生長過程。 Wu 等[5]用Realizedk-ε模型探究了多級流化床的氣固流動特性。 2015 年,Wray 和Agarwal[6]在兩方程的基礎上提出了一方程模型,即Wray - Agarwal(WA)湍流模型,并在單氣相流動模擬中獲得了成功。 Han 等[7]在WA 模型的基礎上提出了一種新的分離渦模型,并應用于壁面有界分離流動的模擬中,結果表明新提出的模型與實驗結果吻合較好。 相對于單相流,多相流的流型復雜多變、相間相互作用強,數學描述難度更大,且計算負荷高。WA 模型對旋流和復雜壁面條件進行了優化,它在復雜流動條件下具有效率高、 成本低、精度高、收斂性好等優點。 Hu 等[8]在旋轉填充床中引入WA模型描述氣液流動,基于文獻的持液量和數均分子量數據驗證了模型的有效性。 Tang 等[9]、Liang等[10]將WA 模型應用于方形鼓泡塔模擬氣液的運動,模擬結果與實驗數據吻合良好。 Shao 等[11]首次將WA 模型應用于循環流化床立管的氣固流動模擬中,結果表明,WA 模型能夠很好地預測立管上部的顆粒分布,但不能很好地預測立管下部的顆粒分布。 在大部分區域WA 模型對顆粒濃度和質量通量的徑向分布的預測總體上優于常規模型,但在中心部分仍存在預測過高的問題。

本文根據氣固兩相流動特征,將WA 湍流模型應用于鼓泡流化床的模擬研究,并將其與Standardk-ε模型和SSTk-ω模型進行對比分析,驗證WA模型的可靠性及適用性,為WA 模型進一步應用于化學反應器內多相流動提供一定基礎。

1 數值模擬方法

1.1 幾何模型及網格劃分

基于Dubrawski 等[12]的實驗,本文建立流化床的三維物理模型。 如圖1 所示,反應器分為流化段和擴展段。 流化段直徑0.133 m,高0.96 m;擴展段直徑0.19 m,高0.4 m。 由于真實的流化床入口分布器結構比較復雜,為網格劃分和數值模擬帶來了一定困難,所以本文將流化床底部簡化為一個均勻進氣的圓面,根據模型在ICEM 中進行了反應器的建模與結構網格的劃分。

圖1 幾何模型及網格劃分Fig.1 The geometric model and meshing

1.2 網格無關性驗證

為了獲得與網格無關的計算結果,本節共生成了29 100、43 008、92 880、226 600 和366 000 這5 個數量的網格,并探究了其對床層壓降的影響。圖2 為床層壓降隨網格數量的變化,從圖中可以看出當網格數量達到92 880 時,床層壓降便趨于穩定,綜合考慮計算成本和精度,后續都采用此網格尺寸。

圖2 網格數目對床層壓降的影響Fig.2 The effect of grid number on bed pressure drop

1.3 數學模型

1.3.1 基本控制方程

本文采用基于歐拉-歐拉方法的雙流體模型(two-fluid model,TFM)來進行鼓泡流化床中氣固兩相流動的數值模擬研究,它將氣固兩相都看作連續相,守恒方程如下。

連續性方程

氣相

式中:p為壓力;Sv為非均相反應導致的動量變化源項;g為重力加速度;β為相間動量交換系數;τ為剪切應力張量。

1.3.2 顆粒動理學理論

雙流體模型將顆粒相按照連續相來處理,顆粒動理學理論解決了如何描述固相黏度和壓力的問題,使得雙流體模型封閉。

顆粒相壓力

式中:Θs為顆粒擬溫度;enm為碰撞恢復系數。

顆粒相剪切黏度

式中:μs,col為顆粒相剪切黏度的碰撞項;μs,kin為顆粒相剪切黏度的動力項;μs,fr為顆粒相剪切黏度的摩擦項。

1.3.3 曳力模型

本文采用Gidaspow 曳力模型,它是Wen-Yu 曳力模型和Ergun 方程的結合,能對粗網格條件下的Geldart B 類顆粒進行準確的模擬。 模型方程如下。

當αg＞0.8 時,由Wen-Yu 模型給出

式中:k為湍流脈動動能;ω為比耗散率;Γ為有效擴散系數;D為交叉擴散系數。

3) WA 模型

Wray-Agarwal(WA)模型是基于k-ω模型導出的新的單方程模型,它將k和ω方程合并為一個R(R=k/ω)運輸方程。

式中:v為運動黏度;ρ為流體密度。

式中:σR、σkε、σkω為湍流普朗特數。

湍流黏度由下式計算:

WA 模型運輸方程中的模型常數取值如下:C1kω=0.082 9,C1kε=0.112 7,σkω=0.72,σkε=1.0,κ=0.41,Cω=8.54,Cm=8.0。

1.4 邊界條件及求解方法

氣相設置為空氣,采用速度入口,設置整個圓柱底面均勻進氣,速度0.4 m/s;出口為壓力出口;壁面為無滑移壁面。 第二相固相設置為顆粒,其密度2 330 kg/m3,均一粒徑103 μm。 固相的初始床層高度為0.8 m,體積分數為0.55。

本文所有的模擬計算均在ANSYS-Fluent 上進行,求解器設置為三維雙精度,并采用多核心節點并行計算。 WA 湍流模型通過軟件提供的接口,以用戶自定義標量方程(User Defined Scalar,UDS)的形式添加到Fluent 中,其他模型參數等通過用戶自定義函數(User Defined Function,UDF)的形式編譯到軟件中進行求解。 入口處R方程的值設置為1 ×10-5,壁面處R的值為0,出口處為0 通量。 同時為了保證計算過程的穩定收斂,R的松弛因子為0.2。求解時采用壓力-速度耦合的Phase Coupled SIMPLE 算法,求解的時間步長設定為0.001 s,每個時間步的最大迭代數為40。 為了保證模擬結果穩定性,總共計算30 s,并取10 ～30 s 內的平均值來進行后續數據的分析。

2 結果與討論

2.1 固含率分布瞬時云圖

圖3 為不同湍流模型模擬得到的床層固含率瞬時分布云圖。 3 種湍流模型得到的床層的最大膨脹高度基本相同,整個床層的近壁面區域固相體積分數較大,這也是由于床層中的大氣泡主要在流化床的中心處產生,使得近壁面的顆粒相對較多。 比較不同時刻下的切面固含率的瞬時分布云圖可知,隨時間的變化,流化床內呈現周期性大氣泡的生長與破裂,符合實際工況,3 種湍流模型均可對氣固流化床的流動進行模擬。

圖3 3 種湍流模型的床層顆粒固含率瞬時云圖Fig.3 Instantaneous cloud of bed particle solid holdup with three turbulence models

2.2 床層壓降分布

在無內構件的空塔情況下,流化床的床層壓降可由下式計算得到

式中:Δp為床層壓降;H0為初始床層高度;ε0為初始床層空隙率。

圖4 為不同湍流模型預測的床層壓降隨時間的變化規律。 由圖可知,流化床在2 s 左右可達到穩定流化狀態,壓降隨大氣泡不斷地生長與破裂呈現有規律的波動。 Standardk-ε和SSTk-ω湍流模型的模擬結果相差不大,絕大部分壓降超過理論值;WA 模型則在理論值附近波動,其預測結果更準確。

圖4 3 種湍流模型預測的床層壓降隨時間的變化Fig.4 Variation of bed pressure drop predicted by three turbulence models with time

表1 為不同湍流模型預測的平均床層壓降模擬值與理論值的對比。 由式(26)計算得到的壓降理論值為5 499.36 Pa,Standardk-ε湍流模型、SSTk-ω湍流模型和WA 湍流模型的預測值的相對誤差分別為6.4%、8.6% 和1.6%。 結果表明,WA模型的預測值的相對誤差更小,模擬結果更準確。

表1 3 種湍流模型預測的床層壓降與理論值的對比Table 1 Comparison between simulated and theoretical values of bed pressure drop predicted by three turbulence models

表2 為采用不同湍動模型反應器內的平均湍動能。 由表可知,通過WA 模型計算的湍動能低于其他兩種湍流模型的模擬值,導致床層壓降模擬值稍小于其他兩種湍流模型的計算結果。

表2 不同湍流模型的整體平均湍動能Table 2 Overall average turbulent kinetic energy of different turbulence models

2.3 軸向空隙率分布

圖5 為3 種不同湍流模型模擬得到的軸向平均空隙率模擬值與實驗值的對比圖。 從圖中可以看出,不同湍流模型的空隙率沿軸向高度的模擬結果規律相近。 SSTk-ω湍流模型的模擬值在流化床底部顯示出與實驗值更好的一致性,而WA 湍流模型在流化床的中部(0.75 m)附近表現出更好的一致性。 總體來說,3 種湍流模型均與實驗結果吻合較好,均可對流化床的氣固流動進行準確地模擬。

圖5 軸向平均空隙率模擬值與實驗值的對比Fig.5 Comparison between simulated and experimental values of axial average voidage

2.4 徑向固含率分布

圖6 為不同湍流模型在4 個軸向高度處的徑向固含率分布模擬值與實驗值的對比。 從圖中可以看出,3 種湍流模型對于徑向固含率的模擬結果的整體走勢是相同的,都呈現一種中間低、兩邊高的對稱分布特征,這與實驗結果相同,也與2.1 節的固含率瞬時分布云圖一致。

圖6 3 種湍流模型在不同高度的徑向固含率分布模擬值與實驗值的對比圖(ug =0.4 m/s)Fig.6 Comparison between simulated and experimental values of radial solid holdup distribution of three turbulence models at different heights (ug =0.4 m/s)

在流化床的底部區域(H=0.24 m),3 種湍流模型都低估了中心處的固相體積分數,而在壁面附近,WA 湍流模型與實驗值更為接近。 在流化床的中部區域(H=0.40 m 及H=0.56 m),WA 湍流模型的預測值與實驗結果在壁面處及流化床中心處均吻合較好。 在流化床的上部區域(H=0.72 m),3 種湍流模型對于徑向固含率的預測趨勢一致,在流化床中心處與實驗值吻合較好,而在壁面附近都高估了徑向固含率。

2.5 顆粒速度分布

顆粒速度在氣固流化床流體力學特征的測定中起著重要的作用,本小節給出一些WA 模型的模擬結果。 圖7 為顆粒軸向速度時均分布云圖,圖8 為不同床層高度處顆粒軸向速度的徑向分布。 從圖中可以看出,顆粒整體上呈現一種中間區域向上運動、近壁面區域向下運動的趨勢。 這種流動規律與氣泡的形成、運動和分裂有關。 由于邊壁效應的存在,中心區域的阻力較小,帶著顆粒的氣泡在流化床的中心向上移動,移動到一定高度時,氣泡破裂并釋放出其中的顆粒,然后粒子在重力的作用下沿著壁面向下運動,從而形成穩定的流化狀態,也導致了內部顆粒的非均勻徑向分布。

圖7 顆粒軸向速度時均分布云圖Fig.7 Time-averaged distribution cloud of particle axial velocity

圖8 不同床層高度處顆粒軸向速度的徑向分布Fig.8 Radial distribution of particle axial velocity at different bed heights

3 結論

本文成功地將WA 湍流模型應用于氣固流化床流動的模擬中,并與經典的兩方程模型進行對比,驗證了WA 模型的可靠性與適用性。 結果表明:

(1)WA 湍流模型對床層壓降的預測值與理論值的相對誤差最小,僅為1.6%。

(2)通過對平均軸向空隙率的分析可知,3 種湍流模型的預測值均在誤差允許的范圍內,SSTk-ω湍流模型在床層底部與實驗結果有更好的一致性,WA 湍流模型在床層的中部附近吻合得更好。

(3)通過對不同高度處平均徑向固含率分布的分析可知,在流化床的中部區域,即H=0.40 m、H=0.56 m 附近,WA 湍流模型的預測值在整體上與實驗值吻合最佳;在流化床的上部區域WA 湍流模型的預測值也比Standardk-ε和SSTk-ω湍流模型的預測值更接近實驗值。