LBM-DEM四向耦合噴動(dòng)床內(nèi)介尺度模擬與分析

， ， ， ， ，

(華北電力大學(xué) 動(dòng)力工程系，保定 071003)

1 引 言

噴動(dòng)床具有氣固接觸時(shí)傳熱效率高、流形簡(jiǎn)單、操作壓降低和便于處理熱敏物質(zhì)等優(yōu)點(diǎn)[1]。研究噴動(dòng)床內(nèi)氣固兩相流，有助于揭示噴動(dòng)床內(nèi)顆粒劇烈而復(fù)雜的運(yùn)動(dòng)機(jī)理。研究?jī)上嗔鞯姆椒ㄖ饕袑?shí)驗(yàn)方法、理論分析和數(shù)值模擬。實(shí)驗(yàn)方法往往對(duì)實(shí)驗(yàn)設(shè)備要求很高，實(shí)驗(yàn)成本也較大，而且不能從微觀上揭示顆粒的運(yùn)動(dòng)特性。理論分析一般只應(yīng)用于通過(guò)合理簡(jiǎn)化能得到比較簡(jiǎn)單模型的體系，但隨著體系復(fù)雜度的增加，其建模和求解難度會(huì)很大。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，數(shù)值模擬克服了實(shí)驗(yàn)方法和理論分析的缺點(diǎn)，以其成本低和效率高的特點(diǎn)，引起了眾多學(xué)者的關(guān)注。目前，用于模擬氣固兩相流的方法主要有歐拉-歐拉法和歐拉-拉格朗日法。其中，氣相的求解都需要根據(jù)有限容積法等求解復(fù)雜的N-S方程，邊界難以處理，且不易于并行計(jì)算。近些年發(fā)展起來(lái)的格子Boltzmann方法能夠有效地克服上述傳統(tǒng)數(shù)值模擬方法的缺陷[2]。

格子Boltzmann方法是20世紀(jì)80年代中期提出的一種從介觀角度出發(fā)并基于分子動(dòng)理論的數(shù)值模擬方法[3]。Filippova等[4]建立了格子Boltzmann方法與離散單元法的耦合模型，但采用單向耦合，只考慮了流體對(duì)顆粒的作用，忽略了顆粒對(duì)流體的作用。而Cui等[5,6]也主要分析了流體對(duì)顆粒的影響，顆粒對(duì)流體的作用分析較少。Govan等[7]雖然考慮了顆粒對(duì)流體的作用，但卻只采用了雙向耦合。Kruggel-Emden等[8]僅研究了顆粒和顆粒群對(duì)流體的作用。Lin[9]用三種不同的模型研究了鼓泡床中的氣泡行為，沒(méi)有深入探討顆粒相的運(yùn)動(dòng)。Wang等[10]提出將格子Boltzmann方法與離散單元法耦合的離散顆粒模型，但其所采用的硬球模型將顆粒之間的碰撞視為瞬間剛性碰撞，無(wú)法給出碰撞的細(xì)節(jié)，且處理兩個(gè)以上顆粒之間的碰撞問(wèn)題具有一定的局限性。

本文將格子Boltzmann方法與離散單元法相結(jié)合，通過(guò)建立的介觀尺度氣固兩相流LBM-DEM四向耦合模型，采用Fortran語(yǔ)言編程，以研究顆粒在噴動(dòng)床內(nèi)的流動(dòng)特性作用機(jī)理，為噴動(dòng)床的設(shè)計(jì)和優(yōu)化提供一定的參考。

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

基于LBM-DEM四向耦合的數(shù)學(xué)模型，考慮到固體運(yùn)動(dòng)對(duì)流場(chǎng)的影響，并參照文獻(xiàn)[10]的處理方法，采用修正的格子Boltzmann方法模擬流場(chǎng)，采用Boltzmann理論計(jì)算控制體內(nèi)流體受顆粒的作用力，將時(shí)驅(qū)硬球模型改為離散單元法軟球模型處理顆粒間以及顆粒與壁面間的碰撞問(wèn)題，將EMMS模型[11]改為Gidaspow模型[12]來(lái)計(jì)算顆粒所受曳力。

2.1 氣相運(yùn)動(dòng)模型

標(biāo)準(zhǔn)的格子Boltzmann方法一般用來(lái)求解單向流體的流動(dòng)。而本文考慮顆粒運(yùn)動(dòng)對(duì)流體流動(dòng)的影響，所以需要在格子Boltzmann方程中加入附加碰撞相：

(1)

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

(2)

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

(3)

式中εs為顆粒體積分?jǐn)?shù),τ為無(wú)量綱松弛時(shí)間。

(4)

式中εg為氣體的空隙率。

將大渦模擬(Large Eddy Simulation)應(yīng)用到LBM可以模擬高雷諾數(shù)的流動(dòng)[14]，即用湍流粘度系數(shù)υe來(lái)反映亞網(wǎng)格的小尺度渦對(duì)流動(dòng)結(jié)構(gòu)的影響。在Smagorinsky亞網(wǎng)格尺度應(yīng)力模型[15]中，υe由應(yīng)變速率決定,

(5)

式中Cs為Smagorinsky常數(shù),Si j為應(yīng)變速率張量。由文獻(xiàn)[16]得,

(6)

式中ek i為離散速度ei在k方向上的分量,Qi j為動(dòng)量通量。

通過(guò)τt=τ+τe代替式(2)的τ來(lái)實(shí)現(xiàn)LES應(yīng)用到LBM，其中τe為有效碰撞渦松弛時(shí)間，則

(7)

由式(6，7)得到一個(gè)一元二次方程，求解得到修正后的松弛時(shí)間τt為

(8)

2.2 顆粒運(yùn)動(dòng)模型

采用牛頓第二定律描述顆粒的運(yùn)動(dòng)，單顆粒i的受力運(yùn)動(dòng)可以表示為

(9)

式中Fy,i為顆粒所受曳力，F(xiàn)P為顆粒所受周圍顆粒的碰撞力。

顆粒之間的碰撞力由彈性力和阻尼力組成。法向力fn,ij分解為法向彈性力fcn,ij和法向阻尼力fdn,ij，切向力ft,ij分解為切向彈性力fct,ij和切向阻尼力fdt,ij[17]。

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

顆粒所受法向力：

fcn,ij(t)=fcn,ij(t-Δt)-knΔδn

(10)

(11)

顆粒所受切向力：

fct,ij(t)=fct,ij(t-Δt)-ktΔδt

(12)

(13)

(14)

2.3 氣固之間的耦合

采用牛頓第三定律解決氣固之間的耦合問(wèn)題。在一個(gè)控制體內(nèi)，采用Boltzmann理論求解顆粒對(duì)流體的作用力[18]，用Gidaspow模型求解流體對(duì)顆粒的作用力。

顆粒對(duì)流體的作用力可表示為

(15)

流體對(duì)顆粒的作用力可表示為

(16)

(17)

顆粒雷諾數(shù)

(18)

為了減小二維系統(tǒng)與三維系統(tǒng)的差異，本文采用文獻(xiàn)[19,20]提出的方法來(lái)計(jì)算空隙率，

(19)

式中 二維系統(tǒng)空隙率

(20)

控制體內(nèi)流固之間的相互作用力以Boltzmann理論計(jì)算出的力為基準(zhǔn)，因此要對(duì)Gidaspow模型計(jì)算出的曳力進(jìn)行修正，來(lái)滿足牛頓第三定律。計(jì)算標(biāo)度因子為

(21)

通過(guò)標(biāo)度因子對(duì)控制體內(nèi)每一個(gè)顆粒受到的曳力進(jìn)行修正：

(22)

利用上述模型，采用文獻(xiàn)[21]的模擬參數(shù)，首先模擬了一段時(shí)間內(nèi)射流在噴動(dòng)床中的演化過(guò)程。整個(gè)過(guò)程與文獻(xiàn)[21]的結(jié)果及文獻(xiàn)[10]觀察到的現(xiàn)象吻合較好。

3 模擬對(duì)象及參數(shù)

模擬對(duì)象在無(wú)量綱下為25×200×0.5的單孔射流準(zhǔn)三維矩形噴動(dòng)床，噴口寬度在無(wú)量綱下為2，模擬對(duì)象如圖1所示。顆粒相采用2970個(gè)直徑為135 μm的球型顆粒，密度為2500 kg·m-3。氣相密度為1.1795 kg·m-3，氣體粘度為1.8872×10-5kg·(m·s)-1。具體模擬參數(shù)列入表1。

本文模擬二維矩形噴動(dòng)床內(nèi)顆粒在底部中間射流作用下的運(yùn)動(dòng)過(guò)程。氣體在壁面處設(shè)置為無(wú)滑移邊界條件，顆粒在壁面處設(shè)置為滑移邊界條件。

圖1 模擬對(duì)象幾何尺寸(無(wú)量綱)

Fig.1 Geometry size of the simulated object (dimensionless)

4 模擬結(jié)果分析

4.1 模擬序列圖

在噴口速度為1.1 m/s時(shí)，0～1 s內(nèi)顆粒流化過(guò)程如圖2所示。隨著射流從床層底部中心位置的噴口處不斷射入，床中逐漸形成的氣泡運(yùn)動(dòng)到床層頂部并發(fā)生破裂。顆粒在中間主氣流的曳力作用下由床層底部向上運(yùn)動(dòng)，運(yùn)動(dòng)到床層表面，伴隨著氣泡的破裂開(kāi)始向兩邊拋灑。兩側(cè)環(huán)隙區(qū)氣流微弱，顆粒主要受自身的重力作用向床層底部運(yùn)動(dòng)。而床層底部顆粒受到中間氣流的卷吸作用由兩邊向中間匯聚，運(yùn)動(dòng)到中間位置又隨著主氣流重新向上運(yùn)動(dòng)。以此循環(huán)，形成了噴動(dòng)床中的環(huán)核流動(dòng)結(jié)構(gòu)。

表1 模擬參數(shù)

Tab.1 Simulation parameters

類別 模擬參數(shù)有量綱無(wú)量綱氣體床寬/m0.0067525床高/m0.054200進(jìn)氣口寬度/m5.4×10-42氣體密度/(kg·m-3)1.8872×10-51氣體粘度/(kg·(m·s)-1)1.8872×10-55.27×10-3格子長(zhǎng)度/m2.7×10-41時(shí)間步長(zhǎng)/s2.4×10-51顆粒彈性恢復(fù)系數(shù)0.90.9摩擦系數(shù)0.10.1顆粒密度/(kg·m-3)25002119.5顆粒直徑/m1.35×10-40.5顆粒數(shù)目29702970時(shí)間步長(zhǎng)4.8×10-60.2

圖2 模擬結(jié)果隨時(shí)間序列圖

Fig.2 Simulation results with time series

4.2 顆粒速度分布

0～1 s內(nèi)，顆粒在不同射流速度和床層高度下水平方向速度的平均值如圖3所示。顆粒水平方向速度大致呈左右中心對(duì)稱，在床層較低位置的顆粒由于中間射流的卷吸作用由兩側(cè)向中間匯聚；在床層較高位置，伴隨著氣泡的破裂，氣體能量逐漸耗散，顆粒所受曳力逐漸減小，速度也逐漸減小。當(dāng)射流速度增加后，在相同床層高度處的顆粒速度變大，運(yùn)動(dòng)更加劇烈。

0～1 s內(nèi)，顆粒在不同射流速度和床層高度下豎直方向速度的平均值如圖4所示。在中間噴動(dòng)區(qū)域，由于氣流速度較大，顆粒所受曳力大于自身重力，速度向上，且隨著高度的增加，氣流的拖曳作用減弱，顆粒的速度也逐漸減小。而在兩側(cè)區(qū)域，氣流速度較小，顆粒所受曳力小于自身重力，顆粒速度向下。當(dāng)射流速度增加后，中間噴動(dòng)區(qū)和兩側(cè)環(huán)隙區(qū)內(nèi)的顆粒速度均增加。

4.3 氣體速度分布

0～1 s內(nèi)，氣體在不同射流速度和床層高度下水平方向速度的平均值如圖5所示。氣體速度大致呈左右中心對(duì)稱。氣流從床層底部噴口射入，由中間向兩側(cè)發(fā)散，且隨高度的增加，速度逐漸降低。當(dāng)射流速度增加后，在相同床層高度處氣體速度也增加。

0～1 s內(nèi)，氣體在不同射流速度下豎直方向速度的平均值如圖6所示。兩側(cè)環(huán)隙區(qū)氣體速度幾乎為0，對(duì)顆粒作用較小，中間噴動(dòng)區(qū)域氣體速度向上，且速度較大，對(duì)顆粒的卷吸能力較強(qiáng)，同樣隨著床層高度的增加，氣體速度逐漸減小。當(dāng)射流速度增加后，中間噴動(dòng)區(qū)氣體的速度增加明顯，而對(duì)兩側(cè)環(huán)隙區(qū)氣體速度的影響較小，速度幾乎為0。

4.4 床層膨脹高度

在0～1 s內(nèi)，顆粒平均體積分?jǐn)?shù)沿床高的分布如圖7所示。本文認(rèn)為顆粒平均體積分?jǐn)?shù)為0時(shí)對(duì)應(yīng)的床層高度為床層膨脹高度。隨著床層高度的增加，床內(nèi)氣泡逐漸增多，顆粒的體積分?jǐn)?shù)逐漸減小。當(dāng)射流速度由0.9 m/s增加到1.1 m/s時(shí)，顆粒所受曳力變大，運(yùn)動(dòng)更加劇烈，床層高度由0.012 m提高到0.014 m。

4.5 床寬對(duì)顆粒速度的影響

考慮床寬對(duì)噴動(dòng)作用的影響，將床寬分別增加到原床寬的兩倍和三倍進(jìn)行對(duì)比分析。當(dāng)射流速度為1.1 m/s，床層高度為0.0108 m時(shí)，在0～1 s內(nèi)，不同床寬下顆粒水平方向和豎直方向速度的平均值分別如圖8和圖9所示。床寬增加后，顆粒水平方向與豎直方向速度的變化趨勢(shì)幾乎不變，但速度的峰值明顯降低。這是由于床寬增加導(dǎo)致左右壁面對(duì)中間噴動(dòng)區(qū)射流擴(kuò)散的約束能力降低，射流可以更容易向兩側(cè)環(huán)隙區(qū)擴(kuò)散，顆粒所受曳力減小，運(yùn)動(dòng)能力減弱，水平方向與豎直方向速度的峰值降低。

圖3 顆粒水平方向速度

Fig.3 Particle’s horizontal velocity

圖4 顆粒豎直方向速度

Fig.4 Particle’s vertical velocity

圖5 氣體水平方向速度

Fig.5 Gas’s horizontal velocity

圖6 氣體豎直方向速度

Fig.6 Gas’s vertical velocity

圖7 顆粒體積分?jǐn)?shù)沿床高的分布

Fig.7 Distribution of particle volume fraction along bed height

圖8 不同床寬下顆粒水平方向速度

Fig.8 Particle’s horizontal velocity in different bed width

圖9 不同床寬下顆粒豎直方向速度

Fig.9 Particle’s vertical velocity in different bed width

5 結(jié) 論

(1) 氣流從床層底部噴口射入，形成氣泡，并不斷上升，在床層頂部破裂。中間區(qū)域的顆粒向上運(yùn)動(dòng)，到達(dá)床層頂端后向兩邊拋灑，在兩側(cè)環(huán)隙區(qū)顆粒向下運(yùn)動(dòng)，而床層底部顆粒向中間匯聚，從而實(shí)現(xiàn)了顆粒的內(nèi)循環(huán)。

(2) 顆粒水平方向速度隨床層高度的變化波動(dòng)明顯，而豎直方向速度卻呈現(xiàn)出中間高、兩邊低的特點(diǎn)，在中間噴動(dòng)區(qū)，顆粒速度向上，在兩側(cè)環(huán)隙區(qū)，顆粒速度向下；氣體水平方向速度分布呈現(xiàn)中心對(duì)稱，而豎直方向氣體中間區(qū)域速度較大且向上運(yùn)動(dòng)，兩側(cè)區(qū)域速度較小，幾乎為0。

(3) 顆粒平均體積分?jǐn)?shù)沿床高逐漸降低，射流速度增大后，顆粒的運(yùn)動(dòng)更加劇烈，床層膨脹高度增加。

(4) 隨著床寬的增加，顆粒速度的峰值降低。