， ， ， ， ，

(華北電力大學 動力工程系，保定 071003)

1 引 言

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

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

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

2 數(shù)學模型

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

2.1 氣相運動模型

標準的格子Boltzmann方法一般用來求解單向流體的流動。而本文考慮顆粒運動對流體流動的影響，所以需要在格子Boltzmann方程中加入附加碰撞相：

(1)

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

(2)

權函數(shù)B(εs,τ)和結構修正因子ω(εg)的表達式為

(3)

式中εs為顆粒體積分數(shù),τ為無量綱松弛時間。

(4)

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

將大渦模擬(Large Eddy Simulation)應用到LBM可以模擬高雷諾數(shù)的流動[14]，即用湍流粘度系數(shù)υe來反映亞網格的小尺度渦對流動結構的影響。在Smagorinsky亞網格尺度應力模型[15]中，υe由應變速率決定,

(5)

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

(6)

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

通過τt=τ+τe代替式(2)的τ來實現(xiàn)LES應用到LBM，其中τe為有效碰撞渦松弛時間，則

(7)

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

(8)

2.2 顆粒運動模型

采用牛頓第二定律描述顆粒的運動，單顆粒i的受力運動可以表示為

(9)

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

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

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

顆粒所受法向力：

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 氣固之間的耦合

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

顆粒對流體的作用力可表示為

(15)

流體對顆粒的作用力可表示為

(16)

(17)

顆粒雷諾數(shù)

(18)

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

(19)

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

(20)

控制體內流固之間的相互作用力以Boltzmann理論計算出的力為基準，因此要對Gidaspow模型計算出的曳力進行修正，來滿足牛頓第三定律。計算標度因子為

(21)

通過標度因子對控制體內每一個顆粒受到的曳力進行修正：

(22)

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

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

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

本文模擬二維矩形噴動床內顆粒在底部中間射流作用下的運動過程。氣體在壁面處設置為無滑移邊界條件，顆粒在壁面處設置為滑移邊界條件。

圖1 模擬對象幾何尺寸(無量綱)

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

4 模擬結果分析

4.1 模擬序列圖

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

表1 模擬參數(shù)

Tab.1 Simulation parameters

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

圖2 模擬結果隨時間序列圖

Fig.2 Simulation results with time series

4.2 顆粒速度分布

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

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

4.3 氣體速度分布

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

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

4.4 床層膨脹高度

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

4.5 床寬對顆粒速度的影響

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

圖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 顆粒體積分數(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 結 論

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

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

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

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