基于雙流體模型的氣液兩相流三維瞬態(tài)數(shù)值模擬

基于雙流體模型的氣液兩相流三維瞬態(tài)數(shù)值模擬

氣液兩相流模擬鼓泡塔檢測(cè)器具有結(jié)構(gòu)簡(jiǎn)單、操作方便、相距面大、傳熱傳遞質(zhì)量高、傳熱傳遞質(zhì)量好等特點(diǎn)。廣泛應(yīng)用于土木工程、生物工程、環(huán)境和能源工業(yè)領(lǐng)域。過(guò)去20年,有關(guān)鼓泡塔已經(jīng)做了大量的研究工作,包括數(shù)值模擬和實(shí)驗(yàn)測(cè)量,但該類反應(yīng)器的設(shè)計(jì)與放大仍然很困難。由于氣體分散相,鼓泡塔內(nèi)的氣/液相流動(dòng)十分復(fù)雜,到目前為止,文獻(xiàn)報(bào)道的研究結(jié)果大多以宏觀數(shù)據(jù)為主,目前的實(shí)驗(yàn)技術(shù)很難測(cè)量鼓泡塔內(nèi)的流場(chǎng)細(xì)節(jié)。當(dāng)氣含率較高,流型為非均相流時(shí),流場(chǎng)細(xì)節(jié)的測(cè)量則更為困難,因此,高表觀氣速下的局部氣含率及氣液相速度的實(shí)驗(yàn)數(shù)據(jù)報(bào)道很少。但僅憑整體宏觀特性數(shù)據(jù)進(jìn)行鼓泡塔的設(shè)計(jì)與放大很不可靠,整體屬性不能敏銳地反映鼓泡塔的尺寸和結(jié)構(gòu)的影響,只能通過(guò)局部屬性來(lái)表征,因此,有必要獲取鼓泡塔內(nèi)的局部流場(chǎng)信息。如前所述,靠實(shí)驗(yàn)手段測(cè)量局部流場(chǎng)相當(dāng)困難,隨著計(jì)算流體力學(xué)的發(fā)展、計(jì)算機(jī)性能和數(shù)值模型可靠性的提高,依靠數(shù)值模擬獲取更多的局部流場(chǎng)信息已經(jīng)成為可能。本文將采用數(shù)值模擬對(duì)鼓泡塔流動(dòng)進(jìn)行三維瞬態(tài)分析,以獲取更多的流場(chǎng)信息。鼓泡塔數(shù)值模擬屬于氣液兩相流模擬。兩相流模擬主要有兩種方法,即雙流體模型和歐拉-拉格朗日法。其中歐拉-拉格朗日法需跟蹤每個(gè)氣泡,對(duì)計(jì)算資源的要求很高,目前只能模擬低氣含率的兩相流;雙流體模型對(duì)各相分別求解動(dòng)量方程和組分方程,降低了計(jì)算量。本文采用雙流體模型。目前雙歐拉法進(jìn)行氣液兩相流模擬存在兩大困難,一是正確描述兩相間的相互作用力;二是采用恰當(dāng)?shù)耐牧髂Ｐ湍M氣液相湍流。關(guān)于相間作用力的模擬存在很大的爭(zhēng)議,但通常認(rèn)為曳力作用最為主要,Dudukovic等、Krishna等僅考慮曳力就獲得了與實(shí)驗(yàn)相一致的時(shí)均結(jié)果。在曳力的描述中,氣泡尺寸最為關(guān)鍵,然而大部分文獻(xiàn)都采用單一氣泡尺寸模型,這顯然與實(shí)際情況不符。盡管通過(guò)選擇合適的氣泡尺寸可以獲得與實(shí)驗(yàn)結(jié)果相吻合的時(shí)均結(jié)果,但是需要通過(guò)反復(fù)的嘗試、篩選,明顯降低了數(shù)值模擬的預(yù)測(cè)性。為了考慮氣泡的尺寸變化,Ishii等[5,6,7,8,9,10,11,12]已經(jīng)做了一定的工作,采用了不同的氣泡尺寸模型。其中,Ishii等采用氣泡界面?zhèn)鬏斈Ｐ蛠?lái)模擬氣泡界面密度,而Sanyal等采用氣泡簇平衡模型(PBM模型)預(yù)測(cè)氣泡尺寸。氣泡界面?zhèn)鬏斈Ｐ椭恍杞庖唤M氣泡界面密度標(biāo)量方程,而PBM模型需對(duì)每一組氣泡求解組分標(biāo)量方程,顯然,后者比前者需要更多的計(jì)算資源與計(jì)算時(shí)間,特別是PBM模型的各氣泡組分方程相互耦合,增加了收斂的難度。盡管最新的FLUENT商業(yè)軟件中已經(jīng)耦合了該P(yáng)BM模型,但考慮到三維模擬的計(jì)算量,本文仍采用Ishii的單方程氣泡界面密度模型。計(jì)算表明,該模型能大大降低對(duì)計(jì)算量的需求,并具有很好的收斂性,且計(jì)算結(jié)果與實(shí)驗(yàn)值比較吻合。液相湍流的模擬是多相流模擬的難點(diǎn),不同的學(xué)者采用了不同的湍流模型。Ranade等采用單相標(biāo)準(zhǔn)k-ε模型,氣相假定為層流。Sanyal等采用了考慮氣泡影響的k-ε模型模擬液相,氣相湍動(dòng)描述則基于Tchen的理論。Zhou等采用二階矩模型模擬氣液相湍流,并考慮了兩相之間的脈動(dòng)相關(guān)。Zhang等采用大渦模擬對(duì)液相湍流進(jìn)行模擬。可見,并沒(méi)有普遍適用的湍流模型可供使用。鑒于此,本文仍采用修正k-ε模型,氣相湍流模擬則基于Tchen的理論,重點(diǎn)考察單方程氣泡模型能否預(yù)測(cè)鼓泡塔內(nèi)的氣液界面濃度。1數(shù)學(xué)模型1.1湍流黏度模型雙流體模型控制方程可在單相流控制方程的基礎(chǔ)上采用集平均方法推導(dǎo)。本文中假定液體為連續(xù)相,氣體為分散相不可壓縮,且不考慮相間質(zhì)量和熱量傳遞,控制方程如下:質(zhì)量守恒方程?αi?t+??(αiUi)=0(1)?αi?t+??(αiUi)=0(1)動(dòng)量守恒方程?αiρiUi?t+??(αiρiUiUi)=-αi??pi+??(αiτi)+Fi+αiρig(2)?αiρiUi?t+??(αiρiUiUi)=?αi??pi+??(αiτi)+Fi+αiρig(2)τi=(μt+μl)[?Ui+(?Ui)Τ-23Ι??Ui](3)τi=(μt+μl)[?Ui+(?Ui)T?23I??Ui](3)式中Fi為相間作用力,N;τi表示黏性應(yīng)力張量,Pa;μl表示分子黏度,Pa·s;μt為湍流黏度,Pa·s。液相湍流黏度采用k-ε模型來(lái)計(jì)算。模型方程為μt?l=Cμρlk2ε(4)μt?l=Cμρlk2ε(4)?αlρlkl?t+??(αlρlklUl)=??(αlμt?lσk?kl)+αlGk?l-αlρlεl+αlρlΠk?l(5)?αlρlkl?t+??(αlρlklUl)=??(αlμt?lσk?kl)+αlGk?l?αlρlεl+αlρlΠk?l(5)?αlρlεl?t+??(αlρlεlUl)=??(αlμt?lσε?εl)+αlεlkl(Cε1Gk?l-Cε2ρlεl)+αlρlΠε?l(6)?αlρlεl?t+??(αlρlεlUl)=??(αlμt?lσε?εl)+αlεlkl(Cε1Gk?l?Cε2ρlεl)+αlρlΠε?l(6)式中Πk,l、Πε,l表示分散相引起的湍流附加作用項(xiàng);參數(shù)Cε1、Cε2、Cμ、σε及σk為常數(shù),分別取1.44、1.92、0.09、1.0、1.3。1.2氣泡破碎過(guò)程氣泡尺寸是模擬兩相間動(dòng)量傳遞的關(guān)鍵參數(shù),氣泡在鼓泡塔內(nèi)的上升過(guò)程中將經(jīng)歷多次聚并與破碎,尺寸將發(fā)生顯著變化。當(dāng)氣泡所受到的外界應(yīng)力大于表面張力時(shí),氣泡破碎,而當(dāng)兩個(gè)或多個(gè)氣泡相互碰撞且氣泡間液膜小于一定厚度時(shí)將聚并。根據(jù)Ishii等推導(dǎo)的氣泡數(shù)密度傳輸方程,氣泡的界面密度傳輸標(biāo)量方程可表示為?ρa(bǔ)?t+??(ρa(bǔ)Ug)=Sbk-Sco(7)?ρa(bǔ)?t+??(ρa(bǔ)Ug)=Sbk?Sco(7)d=6α/a(8)式中Sbk表示氣泡破碎導(dǎo)致的氣液界面增加,根據(jù)氣泡破碎理論,氣泡破碎的主要機(jī)理是液相湍流渦碰撞,但只有尺寸小于氣泡的湍流渦的碰撞才有效。氣泡破碎速率取決于氣泡與湍流渦的碰撞頻率,并且氣泡Weber數(shù)必須大于某一臨界值氣泡才會(huì)破碎,Wu等依據(jù)氣體分子動(dòng)理論,推導(dǎo)出氣泡碰撞速率表達(dá)式Sbk=118Cbkut(a2α)(1-WecrWe)1/2exp(WecrWe)(9)Sbk=118Cbkut(a2α)(1?WecrWe)1/2exp(WecrWe)(9)We=ρutdσ(10)ut=√2(εd)1/3(11)式中ut表示處于湍流能譜慣性子區(qū)的湍流渦速度。式(7)中Sco表示氣泡聚并所引起的相界面減少。Wu等認(rèn)為氣泡聚并主要因氣泡之間的碰撞產(chǎn)生,碰撞頻率與湍流脈動(dòng)有關(guān),若只考慮氣泡之間的二元碰撞,聚并速率可以表示為Sco=-Cco(uta2)3π×[1α1/3max(α1/3max-α1/3)]{1-exp[-α1/3max(αmax/α)1/3-1]}(12)Sbk與Cco為可調(diào)節(jié)參數(shù),本文中取0.3;Wecr數(shù)根據(jù)文獻(xiàn)取0.18;αmax為最大氣含率,取0.8。式(11)中耗散率可通過(guò)求解k-ε方程獲取。本文僅考慮兩相之間的曳力作用,曳力系數(shù)采用Schiller&Naumann模型FD?g=-0.75αg(1-αg)ρlCDdB|Ug-Ul|(Ug-Ul)CD={24(1+0.15Re0.687)/ReRe≤10000.44Re>1000(13)2速度進(jìn)口邊界本文模擬一圓柱形鼓泡塔,塔徑D=190mm,凈液位高度H=1000mm,分布器采用多孔分布結(jié)構(gòu),孔徑1mm,網(wǎng)格結(jié)構(gòu)如圖1所示。模擬中,鼓泡塔頂部采用壓力出口邊界條件,底部采用速度進(jìn)口邊界,壁面為非滑移邊界。初始時(shí)間步長(zhǎng)Δt=0.001s,待計(jì)算達(dá)到穩(wěn)定后取Δt=0.005s。離散格式時(shí)間導(dǎo)數(shù)采用一階迎風(fēng)格式,空間導(dǎo)數(shù)采用Quick格式,壓力-速度耦合采用phase-coupledsimple方法,壓力松弛因子為0.7,動(dòng)量方程松弛因子取0.3,其余采用默認(rèn)值。氣泡界面密度傳輸模型通過(guò)在Fluent中添加一組標(biāo)量傳輸方程來(lái)實(shí)現(xiàn),標(biāo)量方程的聚并項(xiàng)和破碎項(xiàng)可通過(guò)用戶自定義程序來(lái)實(shí)現(xiàn)。待流動(dòng)充分發(fā)展并達(dá)到擬穩(wěn)態(tài)后,即可對(duì)特定位置上的瞬時(shí)物理量進(jìn)行記錄,并保存在外部文件中,物理量的實(shí)時(shí)記錄可通過(guò)UDF程序來(lái)實(shí)現(xiàn)。3結(jié)果與討論3.1氣含率分布模擬結(jié)果本文模擬了0.02、0.05、0.10、0.12、0.15、0.20、0.25、0.30m·s-18個(gè)表觀氣速,涉及均相流型與非均相流型。圖2描述了表觀氣速0.12m·s-1時(shí)中心軸截面的液相速度流型,箭頭長(zhǎng)度表示速度大小,顏色表示氣含率。由圖可知,本文采用的模型有效地揭示出鼓泡塔內(nèi)氣液兩相流的瞬態(tài)行為,特別是大渦的發(fā)展,這種大尺度渦的發(fā)展呈現(xiàn)明顯擺動(dòng)。整體氣含率是表征鼓泡塔性能的重要參數(shù),也是衡量數(shù)值模擬結(jié)果是否可靠的重要依據(jù)。圖3虛線表示整體氣含率模擬值,實(shí)線為實(shí)驗(yàn)結(jié)果,可以看出模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本吻合。模擬曲線反映了氣含率隨表觀氣速增加而增大的趨勢(shì),但隨著表觀氣速增加,氣含率的增加幅度減小,若將氣含率表示為表觀氣速的函數(shù),則指數(shù)因子應(yīng)小于1。氣含率分布使鼓泡塔在徑向方向形成密度差,推動(dòng)了整體循環(huán)流動(dòng)的形成,因此,預(yù)測(cè)局部氣含率分布十分重要。圖4描述了表觀氣速為0.12m·s-1時(shí)氣含率的徑向分布,由圖可知,隨著高徑比的增大,塔中心氣含率逐步提高,分布曲線斜率增大,這是因?yàn)殡S著高徑比的增大,流型逐漸發(fā)展,多孔分布器形成的小氣泡經(jīng)過(guò)不斷的聚并,形成越來(lái)越多的大氣泡,并逐步向塔心聚集所致。但高徑比為2、3處的氣含率分布差別已很小,說(shuō)明流動(dòng)已發(fā)展充分。另外,圖4對(duì)氣含率分布模擬值與文獻(xiàn)實(shí)驗(yàn)結(jié)果進(jìn)行了比較。盡管文獻(xiàn)取于H/D=2.8處,而模擬結(jié)果取于H/D=3處,但可以認(rèn)為這兩個(gè)位置上的氣含率基本相同?？梢钥闯?模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好,但塔中心模擬值比實(shí)驗(yàn)值稍高,隨著遠(yuǎn)離中心軸,模擬值反而低于實(shí)驗(yàn)值,但誤差都小于10%,進(jìn)一步靠近壁面,模擬結(jié)果與實(shí)驗(yàn)結(jié)果很接近。另外,模擬曲線相對(duì)光滑,而實(shí)驗(yàn)曲線則顯得過(guò)于突兀。圖5為不同表觀氣速下的氣含率徑向分布,可以看出,局部氣含率隨表觀氣速的增大而增大,特別是塔中心,局部氣含率增加明顯。當(dāng)表觀氣速較高時(shí)(>0.15m·s-1),壁面附近的氣含率增加不明顯,甚至有所降低,這是因?yàn)楸诿嫣幍臍夂手饕怯苫亓鳉馀萁M成,當(dāng)表觀氣速增大時(shí),大氣泡增加,小氣泡減少,導(dǎo)致壁面附近的回流氣泡減少,氣含率降低。另外,表觀氣速很低時(shí)(Ug=0.02m·s-1),氣含率分布相當(dāng)均勻,說(shuō)明此時(shí)正處于均相流型,隨著氣速的增大,氣含率分布曲線越來(lái)越陡,說(shuō)明流型逐步由均相流轉(zhuǎn)變?yōu)榉蔷嗔鳌?.2表觀氣速對(duì)液相速度的影響圖6為中心線軸向液相速度Riquarts關(guān)系式預(yù)測(cè)值與模擬結(jié)果的比較。從圖中可以看出,模擬值與關(guān)系式預(yù)測(cè)值吻合較好。當(dāng)氣速低于0.15m·s-1,模擬值比計(jì)算值要高10%左右,而當(dāng)氣速高于0.15m·s-1,模擬結(jié)果稍低于計(jì)算值,但都處于工程問(wèn)題所能接受的誤差范圍內(nèi)。圖7給出了表觀氣速為0.12m·s-1時(shí)的軸向液相速度徑向分布。由圖可以看出,塔中心處軸向液相速度隨高徑比的增大而增大,但H/D>2以后,速度增加值減小,曲線斜率變化也減小,說(shuō)明H/D>2以后流動(dòng)逐漸充分發(fā)展。另外,該圖還對(duì)模擬結(jié)果與文獻(xiàn)的實(shí)驗(yàn)值進(jìn)行了定量比較,可以看出,模擬結(jié)果定性地反映了速度徑向分布規(guī)律,模擬較合理地預(yù)測(cè)了液相速度拐點(diǎn)。但鼓泡塔的中心模擬值比實(shí)驗(yàn)值略低,進(jìn)一步靠近壁面,誤差有所增大。目前,有關(guān)鼓泡塔數(shù)值模擬的文獻(xiàn)中僅有少部分工作對(duì)氣相速度進(jìn)行了預(yù)測(cè),并與實(shí)驗(yàn)結(jié)果進(jìn)行比較,大多僅預(yù)報(bào)液相速度,主要是由于氣相的速度不易測(cè)量,特別是表觀氣速較大時(shí)(>0.1m·s-1),氣相速度測(cè)量難度非常大。目前對(duì)非均勻兩相流氣相速度數(shù)據(jù)報(bào)道的甚少,因此,本文僅給出表觀氣速為0.12m·s-1時(shí)軸向氣相速度模擬結(jié)果(圖8)。圖8可以根據(jù)氣相對(duì)液相的滑移來(lái)解釋,對(duì)于一定尺寸的氣泡,氣/液相間滑移速度Δu與氣泡在靜止液體中的上升速度有相同的量級(jí),若在圖7的液相速度曲線上疊加一滑移速度,則可獲得圖8。從圖中可以看出,隨著高徑比增大,氣相速度逐漸增大,特別是塔中心(r/R<0.6)區(qū)域氣相速度增加明顯,這是因?yàn)闅馀菰谏仙倪^(guò)程中聚并形成了大氣泡,上升速度加快,而在壁面附近氣泡尺寸變化不大,以小氣泡為主(大氣泡向塔中心運(yùn)動(dòng)),因此壁面附近氣相速度沿軸向沒(méi)有大的變化。圖8中氣相軸向速度都大于零,說(shuō)明從長(zhǎng)時(shí)均結(jié)果來(lái)看,氣相在各個(gè)空間位置上均呈向上運(yùn)動(dòng),但并不表示氣泡沒(méi)有回流,事實(shí)上壁面附近空間點(diǎn)軸向氣相速度的時(shí)間監(jiān)控曲線中仍存在負(fù)值,說(shuō)明仍有氣泡被液體卷吸回流。圖9顯示了不同表觀氣速下時(shí)均軸向液相速度徑向分布的模擬結(jié)果。所有的時(shí)均曲線均呈拋物線形,中心點(diǎn)的液相速度隨表觀氣速的提高而增加。從圖中可以看出,液相速度的拐點(diǎn)發(fā)生在0.6<r/R<0.8范圍內(nèi),這與文獻(xiàn)結(jié)果相符,但拐點(diǎn)位置隨表觀氣速增加而越來(lái)越靠近壁面,這是因?yàn)闅馑僭龃髸r(shí),氣泡的擺動(dòng)更為劇烈。另外,液相速度最低點(diǎn)隨表觀氣速變化并不明顯,這可能是因壁面影響造成的。圖10(a)、(b)給出了液相與氣相法向雷諾應(yīng)力分量ˉu′zu′z的徑向分布結(jié)果。從圖可以看出,氣/液相的軸向脈動(dòng)速度分量最大值不是在塔中心,而是在r/R=0.6附近,即液相速度分布的拐點(diǎn)位置,說(shuō)明該處速度變化劇烈,湍動(dòng)強(qiáng)度大。比較圖10(a)、(b)可以看出,氣相的平均脈動(dòng)速度要略大于液相,這與周力行的結(jié)論一致。由于缺乏可用實(shí)驗(yàn)數(shù)據(jù),本文只對(duì)模擬結(jié)果進(jìn)行定性分析。Degaleesan通過(guò)實(shí)驗(yàn)獲得162mm鼓泡塔內(nèi)表觀氣速Ug=0.3m·s-1時(shí)的液相軸向雷諾應(yīng)力分量,其最大值約為3500cm2·s-2,比本文模擬結(jié)果大1倍(忽略塔徑的影響),特別是塔中心,模擬結(jié)果要比實(shí)驗(yàn)值低一個(gè)數(shù)量級(jí)。Degaleesan等還測(cè)得了140mm鼓泡塔內(nèi)表觀氣速Ug=0.12m·s-1時(shí)的液體軸向雷諾應(yīng)力分量,若不考慮塔徑效應(yīng),本文模擬結(jié)果僅為實(shí)驗(yàn)結(jié)果的1/3,但曲線的形狀與實(shí)驗(yàn)曲線比較接近。總體說(shuō)來(lái),模擬的雷諾應(yīng)力結(jié)果與實(shí)驗(yàn)值相差較大。Chen等采用相同湍流模型得到的雷諾應(yīng)力模擬結(jié)果也同樣偏低,而Zhou等采用二階矩代數(shù)應(yīng)力模型則較好地預(yù)測(cè)了氣/液相雷諾應(yīng)力?？梢姳疚乃捎玫耐牧髂Ｐ筒荒苡行У啬M氣液相湍流。由于k-ε湍流模型為雷諾時(shí)均模型,時(shí)均化將使許多瞬態(tài)的高頻流動(dòng)特性被過(guò)濾掉,模擬得到的僅為擬瞬態(tài)結(jié)果,因此,采用擬瞬態(tài)結(jié)果來(lái)計(jì)算雷諾應(yīng)力將會(huì)偏低。盡管本文所采用的湍流模型低估了湍流脈動(dòng)速度,但該模型卻較為合理地預(yù)測(cè)了液相湍流動(dòng)能,圖11給出了表觀氣速0.12m·s-1時(shí)的湍流動(dòng)能模擬值,并與實(shí)驗(yàn)值進(jìn)行了比較,預(yù)測(cè)結(jié)果與實(shí)驗(yàn)值較為吻合。3.3表觀氣速對(duì)模擬結(jié)果的影響體積傳質(zhì)系數(shù)是評(píng)估反應(yīng)器性能最為重要的參數(shù),對(duì)于快反應(yīng),體積傳質(zhì)系數(shù)(Kla)決定了反應(yīng)器的生產(chǎn)能力及生產(chǎn)效率。Kla是傳質(zhì)系數(shù)和相間比表面積的乘積。關(guān)于傳質(zhì)系數(shù)Kl有不同的理論,如雙膜理論、湍流渦理論等,但各理論之間存在相互矛盾的地方,因此,本文將著重考察氣液界面表面積密度與表觀氣速的關(guān)系。關(guān)于界面密度已有報(bào)道,本文的模擬結(jié)果將與文獻(xiàn)進(jìn)行定量對(duì)比(圖12),由圖12可以看出,表觀氣速較低時(shí)(<0.10m·s-1)時(shí),模擬結(jié)果與文獻(xiàn)比較接近,而