曳力模型對噴動床氣固兩相流數(shù)值模擬結(jié)果影響的介尺度分析

張慶來， 李 斌， 王雨萌， 陳宇翔

(華北電力大學(xué) 能源動力與機(jī)械工程學(xué)院，河北 保定 071003)

0 引 言

顆粒流體兩相體系廣泛存在于自然界及化工、能源、醫(yī)藥、環(huán)保等工業(yè)生產(chǎn)領(lǐng)域,如沙塵暴、藥粉干燥、造粒、煤氣化等。相間復(fù)雜的相互作用機(jī)理一直是兩相體系研究的重點(diǎn),由于實(shí)驗(yàn)方法很難展現(xiàn)各種復(fù)雜工況下相間相互作用的細(xì)節(jié),因此對于兩相流體系的研究多采用數(shù)值模擬方法。相關(guān)研究已經(jīng)表明,顆粒相應(yīng)力在兩相流中的作用是次要的,相間曳力起著更重要的作用[1]。因此,在兩相流的數(shù)值研究中,曳力模型的選取極為關(guān)鍵,在很大程度上直接決定著模擬結(jié)果的準(zhǔn)確性。

噴動床內(nèi)的氣固兩相系統(tǒng)作為稠密顆粒流體兩相體系的典型代表,在近20年,眾多學(xué)者已針對其進(jìn)行了相當(dāng)深入的研究[2-10]。其中,以傳統(tǒng)CFD數(shù)值模擬方法開展的研究中,針對于時空多尺度(包括宏觀尺度(cm級)及介觀尺度(μm級))領(lǐng)域,相間曳力的計算主要采用根據(jù)實(shí)驗(yàn)數(shù)據(jù)關(guān)聯(lián)所得到的經(jīng)驗(yàn)或半經(jīng)驗(yàn)?zāi)Ｐ?此類模型的共同特點(diǎn)是在單個顆粒曳力模型的基礎(chǔ)上引人顆粒體積分?jǐn)?shù)函數(shù)以描述周圍顆粒的影響。典型代表如Syamlal-O’Brien模型、Wen-Yu模型、Gidaspow模型。其中又以Gidaspow曳力模型的應(yīng)用最為廣泛。此外,近些年,格子玻爾茲曼方法(LBM)的快速發(fā)展使得兩相流研究領(lǐng)域得到重要補(bǔ)充。而基于LBM方法所開展的兩相流研究主要針對于介觀尺度,曳力模型的選取則主要在圍繞于Gidaspow模型與EMMS模型之間[11-13]。當(dāng)前,雖然有學(xué)者針對于不同曳力模型的應(yīng)用進(jìn)行了一定的研究,但多采用成熟商業(yè)軟件開展,且多針對于鼓泡床內(nèi)氣泡行為或均質(zhì)流研究[14-16],而對于噴動床具有非均質(zhì)流化過程的研究還缺少必要的補(bǔ)充及更深入分析。尤其在近些年研究中,有學(xué)者指出考慮了介尺度非均一結(jié)構(gòu)影響的EMMS模型在非均質(zhì)兩相流研究中具有相當(dāng)?shù)膬?yōu)越性[17,18]。

因此,為更深入且更直觀的對比不同曳力模型的差異。本文取Syamlal-O’Brien、Gidaspow、EMMS三種曳力模型(后文分別以S、G、E表示),基于開發(fā)的介尺度LBM-DEM耦合流動程序,從介觀層面對不同曳力模型在預(yù)測噴動床稠密氣固兩相體系流動特征方面的差異進(jìn)行分析,并通過實(shí)驗(yàn)驗(yàn)證各模型的準(zhǔn)確性。當(dāng)然,有研究指出沒有模型可以適合所有工況條件[19]。但我們還是可以找到在一定范圍內(nèi)最優(yōu)的選擇。本文研究工作旨在探索不同曳力模型在低氣速工況條件下的噴動床稠密氣固兩相流研究中的適用性,并為后續(xù)研究工作提供參考。

1 數(shù)學(xué)模型

基于介尺度LBM-DEM耦合模型,如圖1所示。以修正的格子玻爾茲曼方法求解流體相;以適用于處理多顆粒同時碰撞與接觸過程的離散單元法軟球模型求解固體相;相間動量耦合的求解分別基于Sylmal-O’Brien、Gidaspow、EMMS曳力模型。

圖1 D2Q9模型下LBM-DEM耦合狀態(tài)Fig.1 LBM-DEM coupling state under D2Q9 model

1.1 流體運(yùn)動模型(LBM)

考慮到顆粒本身及其運(yùn)動過程將對氣相速度場演化產(chǎn)生一定的影響,因此在速度場的格子Boltzmann方程中加入附加碰撞項(xiàng)進(jìn)行修正[12]。具體如下:

修正后的速度場的格子玻爾茲曼方程演變?yōu)?/p>

(2)

式中:τ為松弛因子;u為氣相格子速度;ρ為氣體密度;fi(r,t)為t時刻r處i方向的速度分布函數(shù);fieq(ρ,u)為i方向速度場的平衡態(tài)分布函數(shù)。

權(quán)函數(shù)B(εs,τ)和結(jié)構(gòu)修正因子ω(εg)的表達(dá)式為

(3)

(4)

式中:εs和εg分別為控制體內(nèi)顆粒的體積分?jǐn)?shù)和氣體的孔隙率,且有εs+εg=1。

1.2 顆粒運(yùn)動模型(DEM)

顆粒運(yùn)動過程由牛頓第二定律進(jìn)行描述,單顆粒i的受力運(yùn)動方程可表示為

(5)

式中:Fy,i為顆粒受到流體的曳力;Fp,j為顆粒i與顆粒j的碰撞力;kc表示與顆粒i發(fā)生碰撞的顆粒數(shù)目。

DEM軟球模型將顆粒間的碰撞視為一種彈性碰撞,可用彈簧、阻尼器及偶聯(lián)器的模型表示[6]。

顆粒間的受力由彈性力和阻尼力組成。阻尼力可分為切向阻尼力fdt,ij和法向阻尼力fdn,ij,彈性力可分為切向彈性力fct,ij和法向彈性力fcn,ij。

根據(jù)胡克定律,顆粒間的作用力可表示如下。

顆粒所受法向力:

(6)

(7)

(8)

顆粒所受切向力:

(9)

(10)

(11)

式中:△δ為相對位移增量;△t為時間步長;μ為摩擦系數(shù);k為彈性系數(shù);η為阻尼系數(shù)。

1.3 氣固耦合(LBM-DEM)

牛頓第三定律用于處理氣固耦合問題,顆粒相對氣相的作用力基于格子Boltzmann理論求解;氣相對顆粒相的作用力求解則分別基于不同的曳力模型。

顆粒對氣體的作用力表示如下:

(12)

式中:h為網(wǎng)格步長;ei為速度配置;i為速度方向。

氣體對顆粒的作用力計算過程如下所示。

基于牛頓第二定律,跟蹤每一顆粒并計算顆粒雷諾數(shù):

(13)

控制容積中單個顆粒的曳力計算:

(14)

式中:dp為顆粒直徑;Vp為顆粒體積。

曳力系數(shù)β的計算分別基于如下曳力模型:

Gidaspow曳力模型:

βG-i=

(15)

(16)

EMMS曳力模型:

(17)

ω(εg)=

(18)

Syamlal- O’Brien模型:

(19)

(20)

+0.5(A-0.06Rep)

(21)

A=εg4.14

(22)

(23)

式中:Cd為阻力系數(shù);Vrs為顆粒沉降速率;ug,vp分別為氣體及顆粒速度;A,B為調(diào)優(yōu)參數(shù)。

2 模擬及實(shí)驗(yàn)設(shè)置

2.1 模擬對象及參數(shù)

模擬對象為寬度×高度(W×H)為0.006 75 m×0.054 m的二維矩形噴動床,噴口位于床體底部中間位置,噴口寬度設(shè)置為2個格子長度。幾何對象如圖2所示。

圖2 模擬對象幾何尺寸Fig.2 Geometry of simulated object

具體模擬參數(shù)設(shè)置詳見表1。

表1 模擬參數(shù)

模擬過程中邊界處理均采用雙精度非平衡外推格式,速度邊界為無滑移邊界。本文研究中,床寬與粒徑尺寸比滿足L:D=50,因此在一定程度上可以忽略邊壁效應(yīng)影響[20]。

2.2 實(shí)驗(yàn)對象及參數(shù)

實(shí)驗(yàn)床體由玻璃材料制作,實(shí)驗(yàn)顆粒選用直徑為120～150 μm的石英砂,質(zhì)量流量計控制氣速,采用高速攝像機(jī)進(jìn)行動態(tài)拍攝,整體實(shí)驗(yàn)裝置如圖3所示。

圖3 實(shí)驗(yàn)裝置Fig.3 Experimental device

為提高實(shí)驗(yàn)與模擬過程的一致性,采用相同顆粒堆積高度(7 mm)作為初始狀態(tài),實(shí)驗(yàn)參數(shù)設(shè)置詳見表2所示。

表2 實(shí)驗(yàn)和模擬參數(shù)

3 結(jié)果分析

分別以三種曳力模型對0.7 m/s氣流入射速度條件下的噴動床進(jìn)行了數(shù)值模擬。預(yù)測得到了2.4 s時間床內(nèi)稠密氣固兩相相互作用的全過程,并通過實(shí)驗(yàn)對預(yù)測結(jié)果進(jìn)行了驗(yàn)證。實(shí)驗(yàn)及模擬結(jié)果如圖4所示。

圖4 噴動床內(nèi)顆粒流化序列實(shí)驗(yàn)及模擬結(jié)果(ux=0.7 m/s)Fig.4 Experimental and simulation results of particle fluidization sequence in spouting bed (ux =0.7 m/s)

實(shí)驗(yàn)結(jié)果顯示:氣流從底部中心噴口射入床體后,床層內(nèi)會生成氣泡并帶動顆粒向上運(yùn)動,隨著射流的持續(xù)注入,氣泡不斷增大,床層高度隨之增加,氣泡突破床層后炸裂,顆粒向兩側(cè)拋灑,呈噴泉狀。此過程中,床體兩側(cè)區(qū)域顆粒不斷向底部噴口處匯聚,又在持續(xù)射流的帶動下向上運(yùn)動,從而形成完整的內(nèi)循環(huán)。床內(nèi)整體流化過程可以分成兩個階段,以首個氣泡炸裂時刻為分界點(diǎn),前階段為氣泡演化過程,床內(nèi)未形成完整流化態(tài),顆粒體積分?jǐn)?shù)接近,可視為均質(zhì)流態(tài);后階段為完全流化過程,流化狀態(tài)成熟,噴泉區(qū)及中心射流區(qū)與兩側(cè)區(qū)域間的流化態(tài)差異較大,因此視為非均質(zhì)流態(tài)。

由模擬結(jié)果可知:三種曳力模型都能完整展現(xiàn)噴動床內(nèi)氣泡演化及顆粒流化的全過程,且完成首個氣泡炸裂的時間均在0.12 s左右,即不同曳力模型對氣泡演化的時間影響較小。三種曳力模型下,氣泡初始階段的演化過程趨于一致,但隨著氣泡接近于床層頂部,氣泡形態(tài)的差異性增大。圖中明顯可以看出,S模型下高床層處的氣泡形態(tài)發(fā)生了擾動,與實(shí)驗(yàn)結(jié)果存在較大的偏差,而G模型與E模型下的模擬結(jié)果相近且與實(shí)驗(yàn)結(jié)果吻合度較高。對于完全流化過程,由床內(nèi)顆粒典型流化態(tài)對比可知,E模型在非均勻流化態(tài)的細(xì)節(jié)展示與實(shí)驗(yàn)最為接近,G模型次之,S模型最差。同樣,模擬所得2.4 s完整流化結(jié)果顯示,S模型所呈現(xiàn)的流化過程混亂度最高,擾動較為劇烈,甚至出現(xiàn)了顆粒飛濺現(xiàn)象,而G模型與E模型下的流化過程則相對比較平穩(wěn)。這一方面體現(xiàn)了G、E兩模型在稠密兩相流數(shù)值計算中具有更高穩(wěn)定性,另一方面也表明E模型在展示噴動床內(nèi)相對非均質(zhì)流動過程的細(xì)節(jié)特征中更具優(yōu)勢。

圖5所示為不同曳力模型下首個氣泡炸裂時刻(0.12 s)床內(nèi)氣相、顆粒相速度場矢量分布。

從圖5可以看出,氣相、顆粒相速度場均呈現(xiàn)不同程度的對稱性,相比較而言,G模型下的氣相及顆粒相速度分布最為規(guī)整,且顆粒相速度場在噴泉區(qū)所呈現(xiàn)的典型環(huán)核流動特征也最為清晰。由于本文模擬對象為中心單噴口噴動床,床內(nèi)初始設(shè)置均為理想的中心對稱設(shè)置,因此G模型的結(jié)果與理論結(jié)果相一致,同時與實(shí)驗(yàn)結(jié)果的吻合度也較高。表明G模型對于均質(zhì)流態(tài)的氣泡演化過程的描述具有較高的優(yōu)勢。這與文獻(xiàn)[14]的研究結(jié)果也相一致。

圖6為床層膨脹高度隨時間變化曲線?？梢园l(fā)現(xiàn)噴動床內(nèi)的床層膨脹高度整體呈波動狀態(tài),不同模型下的床層高度變化趨勢趨于一致。其中,S模型下的床層膨脹高度最高,且與實(shí)驗(yàn)值最為接近,而G模型與E模型的結(jié)果則相對偏低。實(shí)際上,由于氣相格子與顆粒相尺寸的限制,模擬過程的噴口寬度略小于實(shí)驗(yàn)值,所以理論上三種曳力模型下的預(yù)測結(jié)果應(yīng)均低于實(shí)驗(yàn)值。因此,S模型所預(yù)測得到的床層膨脹高度實(shí)際上整體偏高。

圖6 床層膨脹高度對比(ux=0.7 m/s)Fig.6 Comparison of bed expansion height (ux=0.7 m/s)

圖7、圖8分別表示不同床層高度上氣體及顆粒豎直方向速度沿床寬分布。如圖中所示,噴動床內(nèi)氣相速度沿床寬方向呈中心對稱分布。床體兩側(cè)氣速較低,在中心區(qū)域達(dá)到峰值。在不同的床層高度(0.002 7 m及0.009 5 m),其分布趨勢相一致,但隨著床層高度的增加,中心區(qū)域的氣速峰值減小,而兩側(cè)區(qū)域氣速增大。從圖中可知,三種曳力模型下的氣速分布趨勢大致相同,但結(jié)果存在一定差異。其中G模型預(yù)測的氣速峰值最高,E模型預(yù)測的結(jié)果則介于G模型與S模型之間。而S模型下的結(jié)果則整體向左側(cè)偏移。與之相對應(yīng),在相同的床層高度,顆粒豎直方向速度沿床寬分布與氣相速度分布嚴(yán)格一致,但不同模型的結(jié)果差異并不與氣相速度場一一對應(yīng)。如圖8所示,在不同床層高度上 S模型預(yù)測的顆粒速度峰值均為最高,但在較低的床層高度(0.002 7 m),G模型的結(jié)果最低,而在較高的床層高度(0.009 5 m),E模型所計算的結(jié)果最低。

圖7 氣體豎直方向速度沿床寬分布Fig.7 Distribution of gas vertical velocity along bed width

圖8 顆粒豎直方向速度沿床寬分布Fig.8 Particle vertical velocity distribution along bed width

實(shí)際上,氣相及顆粒相的速度分布是兩相相互作用后的結(jié)果表現(xiàn),在實(shí)時的曳力展現(xiàn)中還有所欠缺。因此,取相同床層高度對不同曳力模型的曳力計算結(jié)果進(jìn)行更直觀的分析。

圖9所示為豎直方向曳力沿床寬變化。從圖中可知,沿床寬方向上三種曳力模型得到的曳力分布均表現(xiàn)出極高的對稱性,且在床體中心區(qū)域,三種模型都展現(xiàn)出了由于顆粒團(tuán)聚所引起的曳力下降的特征。在不同的床層高度,三種模型的曳力計算結(jié)果有較大差異,以E模型為基準(zhǔn),在0.002 7 m床層高度,S模型計算的平均曳力值比E模型高16.24%,G模型的計算結(jié)果則比E模型低9.38%。在0.009 5 m床層高度上,S模型計算的平均曳力值比E模型高12.94%,G模型的計算結(jié)果則比E模型低14.69%。

圖9 豎直方向曳力沿床寬變化Fig.9 Vertical drag along the width of bed

在床體的近壁面區(qū)域,三種模型下的曳力結(jié)果差異性更為明顯。如圖10近壁區(qū)曳力變化所示,沿床高方向上G模型與E模型計算的曳力值相近,其均值分別為2.25×10-10N、 3.94×10-10N。而S模型計算的曳力值沿床高方向先增大后減小,其均值為3.26×10-10N,是E模型的近10倍。因此,對于S模型的應(yīng)用而言,需格外注意因曳力計算偏高而可能引起的邊壁效應(yīng)影響。

圖10 近壁區(qū)曳力沿床高變化Fig.10 Variation of drag force along bed height in near-wall area

4 結(jié) 論

本文基于開發(fā)的介尺度LBM-DEM耦合流動程序,分別采用Sylmal-O’Brien、Gidaspow、EMMS曳力模型對噴動床內(nèi)稠密氣固兩相體系的流化特征進(jìn)行了數(shù)值模擬。主要研究結(jié)論如下:

(1)三種曳力模型都能完整展現(xiàn)噴動床內(nèi)氣泡演化及顆粒流化的全過程,其中EMMS模型下的整體預(yù)測結(jié)果與實(shí)驗(yàn)結(jié)果最相吻合,模型適用性最高。

(2)Gidaspow模型在均質(zhì)流態(tài)的氣泡演化過程中極具優(yōu)勢,能更好的展現(xiàn)均質(zhì)流的典型特征。并且在較低顆粒體積分?jǐn)?shù)區(qū)域其對氣相速度的影響作用最小。

(3)整體而言,EMMS模型與Gidaspow模型的預(yù)測結(jié)果較為接近,但在噴動床完全流化過程的預(yù)測中EMMS模型在細(xì)節(jié)處的展現(xiàn)要優(yōu)于Gidaspow模型。

(4)Sylmal-O’Brien模型計算的曳力值整體偏高, 在近壁面區(qū)域其沿床高方向所計算的曳力均值可達(dá)Gidaspow模型與EMMS模型計算結(jié)果的14.48倍與8.27倍。因此Sylmal-O’Brien模型在應(yīng)用時需格外注意由于曳力計算偏高而可能引起的邊壁效應(yīng)的影響。

(5)在預(yù)測結(jié)果準(zhǔn)確性及數(shù)值計算穩(wěn)定性方面,Sylmal-O’Brien模型均差于Gidaspow模型與EMMS模型。因此在噴動床稠密氣固兩相體系的數(shù)值研究中Sylmal-O’Brien模型的適用性較低。