連續(xù)操作密相流化床顆粒停留時間分布特性模擬放大研究

蘭斌，徐驥，劉志成，王軍武

（1 中國科學院過程工程研究所多相復雜系統(tǒng)國家重點實驗室，北京100190； 2 中國科學院大學化學工程學院，北京100049； 3 中國科學院綠色過程制造創(chuàng)新研究院，北京100190； 4 中國石油化工股份有限公司上海石油化工研究院，上海201208）

引 言

氣固流化床因具有良好的氣固接觸效率、較高的傳質傳熱效率和氣固處理能力而被廣泛應用于石油催化裂化反應、金屬礦物的焙燒與還原、化石燃料的燃燒、生物質的熱解和氣化以及造粒等過程[1-5]。然而，由于氣體和顆粒間的復雜作用使得人們至今難以對流化床內氣固兩相流流場進行精確的分析，到目前為止工業(yè)流化床反應器的設計和放大還主要通過實驗逐級放大的方法。但是逐級放大實驗成本高、周期長，因此隨著計算機的發(fā)展，研究人員開始采用CFD 數(shù)值模擬的方法對流化床的放大效應進行研究[6-11]。

目前研究流化床模擬放大主要采用以下兩類方法：歐拉-歐拉雙流體模型（two-fluid model，TFM）和歐拉-拉格朗日CFD-DEM 模型及其簡化方法。例如Verma 等[6]、Che 等[7]采用TFM 研究了流化床直徑對氣泡和固體運動、環(huán)-核結構的影響。Couto等[8]采用耦合k-ε 湍流模型的TFM 方法研究了實驗室和半工業(yè)生物質氣化反應器的氣化溫度、氧氣含量和生物質類型對反應的影響，發(fā)現(xiàn)反應器尺寸越大，氣體停留時間越長，氣化反應越充分，CO 和H2的含量越高，合成氣熱值越高。Lu 等[9]采用EMMS曳力模型預測了四種不同尺度的甲醇制取烯烴反應器的流動特性，認為化學反應動力學模型通常是通過對微尺度流化床實驗數(shù)據(jù)進行擬合得到的，可能并不適用于較大尺寸的反應器。Zhang 等[10]采用CFD-DEM 方法證明了在保證流化床底部進口氣流均勻的情況下，增加流化床長度對煤和煤矸石顆粒分離程度的影響并不大。此外，Gu 等[11]采用多相流質點網(wǎng)格法（multiphase particle-in-cell，MP-PIC）對循環(huán)流化床（circulating fluidized bed，CFB）鍋爐內的氣固流動和富氧燃燒過程進行了研究，發(fā)現(xiàn)實驗室和中試流化床顆粒速度場比工業(yè)流化床更加均勻，工業(yè)CFB 鍋爐的CO、NO 和SO2排放量低于實驗室和中試鍋爐，熱量輸入的增加可以提高脫硫效率。

顆粒停留時間分布（residence time distribution，RTD）是表征反應器內氣固混合程度的重要參數(shù)，近幾十年來，研究人員分別利用實驗、理論和數(shù)值模擬手段來研究氣固流化床，特別是密相連續(xù)操作流化床內顆粒停留時間分布情況[12-15]。實驗方面，白書培等[12]發(fā)現(xiàn)在氣速、粒徑、床高、顆粒流量等操作條件中，氣速是影響停留時間分布的主要因素，隨著氣速的增加，顆粒在床內的混合加劇，氣速達到一定值時，RTD 曲線出現(xiàn)多峰。但高巍等[13]認為顆粒流率、床高、粒徑分布是影響顆粒RTD 的主要因素，氣速則是次要因素。顆粒進料流率也是影響顆粒RTD 的重要因素，流率越大，顆粒停留時間越短，顆粒流動向平推流靠近[13]。Matheson 等[14]發(fā)現(xiàn)相比粗顆粒，細顆粒在床內的混合程度更強，粗顆粒停留時間分布更趨于平推流[13]。不同粒徑顆粒在流化床中停留時間存在差異，隨著顆粒粒徑的增加，顆粒平均停留時間增長[15]，同一流化床內粗顆粒比細顆粒停留時間更長[16]。Yagi 等[17]發(fā)現(xiàn)，如果顆粒轉化受化學反應的控制（最有可能發(fā)生在細顆粒轉化中），則完全轉化時間與顆粒粒徑呈正比；如果顆粒轉化受內擴散控制（最有可能發(fā)生在粗顆粒轉化中），則完全轉化時間與顆粒粒徑的平方呈正比。郝志剛等[18]采用多級橫向內構件研究了粒徑分布較寬的顆粒平均停留時間，發(fā)現(xiàn)700 μm 顆粒平均停留時間是108 μm顆粒的5倍。數(shù)值模擬方面，Geng等[19]、Zou 等[20-21]、Hua 等[22]采 用TFM 研 究 了 固 體 流率、氣速、床高、示蹤劑注入時間、取樣頻率等操作條件對顆粒停留時間分布的影響。Zou 等[23]、Hua等[24]發(fā)現(xiàn)相比均勻結構曳力模型，非均勻結構曳力模型預測的顆粒停留時間分布和實驗結果吻合得更好。B?rner 等[25]提出一種基于相似性模型的放大方法，模擬了氣、液、固三相頂部噴霧造粒機內顆粒的停留時間分布，發(fā)現(xiàn)顆粒所穿過的噴霧區(qū)長度和顆粒速度的變化導致顆粒停留時間不同。Zhao等[26]、Lu等[27]采用CFD-DEM 方法模擬了循環(huán)流化床提升管內顆粒的返混行為。Lan 等[28-29]提出了局部返混指數(shù)的概念，采用CPFD 方法研究了循環(huán)流化床中不同流型顆粒返混特性以及提升管出口結構對顆粒停留時間的影響，指出相對于簡單出口和C型出口，含有L 型和T 型出口的提升管內返混程度更為劇烈，停留時間也更長。系統(tǒng)的流化床顆粒停留時間研究進展見文獻[30]。

顆粒形狀[31-35]和多分散性[36-38]對顆粒的流體力學特性和停留時間分布具有重要影響，為了探究多分散流化床的放大規(guī)律，實現(xiàn)流化床的合理放大，本文以非球形、多分散顆粒曳力模型為基礎[39]，采用基于GPU 大規(guī)模并行的粗粒化CFD-DEM 方法[40]，研究連續(xù)進出料流化床顆粒停留時間和流動特性隨流化床尺寸的變化關系。

1 數(shù)學模型

1.1 粗粒化CFD-DEM 方法

本文所采用的粗粒化CFD-DEM 方法是基于Lu等[40]提出的基于能量最小多尺度離散顆粒方法（energy minimization multi-scale-discrete particle method，EMMS-DPM），但與EMMS-DPM 方法有以下幾個區(qū)別。

（1）假設粗顆粒（coarse grained particle,CGP）的密度等于真實顆粒的密度，即粗顆粒質量mCGP=k3mp，其中，mp是真實顆粒質量，k 是粗粒化率，k=dCGP/dp。而原始EMMS-DPM 方法中的粗顆粒內部存在空隙，其空隙率等于最小流化條件下的空隙率。

（2）實驗采用不規(guī)則形狀的石英砂顆粒，對于非球形顆粒與顆粒/壁面之間接觸力的計算，可以通過組合球元構建非球形顆粒幾何模型的方法實現(xiàn)[41]，但需要計算大量球元之間的碰撞，計算量明顯增大。因此，為提高計算速度，在計算顆粒與顆粒/壁面之間相互作用時，仍將顆粒視為球形。顆粒之間碰撞作用的計算采用Peng等[42]提出的方法。

（3）在CFD-DEM 模型中，氣固相間作用力占主導地位[43]。有研究指出[33]，當顆粒球形度Ψ=1 時，原始的Ergun 曳力系數(shù)關聯(lián)式（黏性項系數(shù)為150，慣性項系數(shù)為1.75）可以很好地預測含球形顆粒的流化床的壓降；但對于非球形顆粒（Ψ≠1），只有當顆粒形狀接近球形時，Ergun 關系式才能較為準確地預測其曳力，因此采用Ganser非球形顆粒曳力模型[33,44]來計算單顆粒曳力系數(shù)CD，并將其引入Di Felice 曳力模型中[45]，具體曳力系數(shù)的計算見1.4節(jié)。

（4）在氣固曳力模型中考慮顆粒的多分散性，參照Beetstra 等[46]、Zhou 等[47]對多分散系統(tǒng)曳力的處理方式，對曳力模型進行修正，計算不同粒徑顆粒所受曳力。

由此可以看出在計算顆粒之間碰撞作用時將顆粒視為球形，而在計算氣固相間曳力時模擬真實顆粒，即引入顆粒球形度這一物理量[通過式（30）、式（31）估算得到球形度的值，此值也符合石英砂顆粒球形度的取值范圍]。這是因為本課題組前期的嘗試表明，如果曳力中不引入球形度的影響，將得到完全錯誤的模擬結果，但是如果在顆粒碰撞時考慮球形度，則計算量將大幅度增加。因此，本文采用的是一種兼顧計算準確性和計算效率的折中近似方案。

1.2 控制方程

氣相連續(xù)性方程

式中，ug為氣相速度，m/s;εg為空隙率。

氣相動量方程

其中，氣相應力張量為

式中，NCGP,cell為網(wǎng)格中粗顆粒數(shù)目；Vcell為網(wǎng)格體積，m3;I為單位張量。

顆粒平動方程

式中，uCGP,i為粗顆粒速度，m/s；右側四項分別為粗顆粒所受壓力梯度力、重力、接觸力和曳力，N 為某一時刻與顆粒i同時相互作用的顆粒總數(shù)。

顆粒轉動方程

1.3 顆粒碰撞模型

粗顆粒之間碰撞力的計算采用彈簧-阻尼模型，顆粒所受接觸力等于法向接觸力和切向接觸力之和，即

顆粒i與j的法向接觸力為

式中，νi、νj為顆粒i和顆粒j的泊松比。

式中，Ri、Rj為顆粒i和顆粒j的半徑；ri、rj代表顆粒i和顆粒j的位置矢量。

法向阻尼系數(shù)

粗顆粒恢復系數(shù)

式中，ep為真實顆粒恢復系數(shù)。

顆粒i與顆粒j的相對速度

顆粒i與顆粒j的法向相對速度

顆粒i與顆粒j的切向接觸力

在模擬中還需考慮滑動摩擦力的限制

式中，μf為滑動摩擦系數(shù)。

1.4 氣固曳力模型

本文采用的曳力系數(shù)以Di Felice 曳力系數(shù)為基礎[45]

式中，dave為網(wǎng)格中真實顆粒平均直徑，m;uave為網(wǎng)格中顆粒平均速度，m/s；εp為網(wǎng)格中顆粒體積分數(shù)。

式中，k 為粗粒化率；mi為粗顆粒質量；ui為粗顆粒速度。

其中，單顆粒曳力系數(shù)[33,44]

顆粒Reynolds數(shù)

Stokes形狀因子

式中，ψ 為顆粒球形度；D 為流化床水力直徑，m；dA為等投影面積球當量直徑，m；dV為等體積球當量直徑，m。由于dA、dV很難通過實驗測得，故假設dA=dV=dsauter，dsauter為顆粒Sauter平均直徑。

牛頓形狀因子

在考慮顆粒非球形影響的Di Felice-Ganser 曳力基礎上，引入Beetstra 等[46]的方法考慮顆粒多分散性的影響，則系統(tǒng)中組分j的曳力系數(shù)可表示為

式中，εj代表網(wǎng)格中組分j 的體積分數(shù)；dj代表組分j的顆粒直徑。

因此，粗顆粒所受曳力為

2 流化床放大方式與模擬設置

圖1 三維連續(xù)進出料流化床裝置模擬Fig.1 Simulation schematic diagram of 3D continuously operated fluidized beds

傳統(tǒng)的流化床放大理論大多基于流體力學相似性[48-51]，一般需要通過直接法[52]（測量氣泡特性參數(shù)、最小流化速度、床層膨脹率以及高速攝像、視頻分析、電容和光學探測等）或間接法[53-55]（壓力波動測量）兩種實驗手段對不同尺寸流化床流體力學相似性進行驗證。為此，有些學者提出了包含多種無量綱參數(shù)的相似性放大規(guī)則，其中最為經(jīng)典的是Glicksman 等[56-57]和Horio 等[58]的放大規(guī)則。但這些規(guī)則一般基于一些假設，如忽略顆粒間相互作用力以及化學反應和傳熱效應，同時要求顆粒大小、最小流化速度、氣速等物性隨流化床尺寸同步改變，并不適用于單獨改變流化床尺寸的放大方式。本文所模擬的四個不同的三維連續(xù)進出料流化床如圖1 所示，長度Lb分別為0.07、0.15、0.31、0.63 m，寬度和高度分別為0.06、0.5 m。

為驗證模型的可靠性，采用趙虎[59]得到的長為0.15 m流化床的實驗結果進行驗證。實驗所測顆粒粒度分布如圖2 所示。Volk 等[60]發(fā)現(xiàn)，從破碎氣泡、減少氣體短路和增強氣固接觸效率的角度來講，寬篩分顆粒要比單一粒徑顆粒更好。因此，為獲取不同粒徑顆粒停留時間信息，同時加快計算速度，在實際CFD-DEM 模擬中將此連續(xù)分布離散為三種粒徑。首先在連續(xù)粒徑分布上取出有限個離散點，求出低階矩[61]，然后使用P-D（product-difference）算法求解得到代表粒徑及其權重（質量分數(shù)）[62]，離散結果如表1所示。

表2 總結了模擬參數(shù)，其中顆粒球形度是根據(jù)實驗測得的最小流化速度結合Hua等[33]、Wen等[63]的研究，由式（30）、式（31）估算得到，DEM 模型中涉及的碰撞參數(shù)均設置為經(jīng)驗值。為獲得足夠的尺度信息，網(wǎng)格尺寸要足夠精細[64-66]，因此選用規(guī)則的六面體網(wǎng)格，尺寸設置為5 mm，約為最大真實顆粒直徑的8 倍。對于Geldart B 類顆粒，該網(wǎng)格大小足以獲得與網(wǎng)格無關的模擬結果[67]。

圖2 實驗所用顆粒尺寸分布Fig.2 Particle size distribution used in experiment

表1 模擬所用顆粒代表粒徑及其質量分數(shù)Table 1 Representative particle sizes and their mass fractions used in simulations

3 結果與討論

3.1 床層壓力分布

不同尺寸流化床壓力時均值軸向分布曲線如圖3所示。可以看出0.15 m 流化床模擬結果和實驗數(shù)據(jù)吻合較好，床內某一高度的壓力隨著流化床長度的增加變化不大，特別是0.31 m 與0.63 m 流化床床層壓降曲線幾乎重合，這是因為當流化床高度不變時，長度的增加對壓降的影響很小。

表2 模擬參數(shù)和氣體、顆粒物性Table 2 Simulation parameters for fluidized bed and properties of gas and particles

圖3 不同長度流化床壓力軸向分布Fig.3 Axial pressure distribution of fluidized beds with different lengths

表3為預測的不同大小流化床的平均存料量和全床時均壓降，從表中可以看出，0.15 m流化床全床壓降模擬結果與實驗結果相對誤差約為6.7%，說明所提出的非球形顆粒曳力模型能夠較為準確地預測床內流動特性。另外系統(tǒng)流動穩(wěn)定后平均存料量隨流化床尺寸的增加而增大，并呈線性關系。全床壓降隨流化床長度的增加變化不大，這與圖3 軸向壓力分布的結果一致。

表3 不同尺寸流化床全床壓降Table 3 Total pressure drop of fluidized beds with different sizes

3.2 固相濃度分布

圖4 為1000 s 時不同尺寸流化床中心豎直截面處固相濃度分布情況。可以看出隨著流化床長度的增加，氣泡數(shù)量增多，直徑變大。0.15、0.31、0.63 m 流化床床層膨脹高度基本一致，而最小床床層膨脹最高，這是由于壁面效應在較小尺寸流化床中更為明顯[68]。床徑較小時，床內氣泡以節(jié)涌的形式向上移動，形成較大的氣栓；而當床徑較大時，氣泡由于聚并而增大，上升速度加快，到達床層表面時發(fā)生破碎，氣泡停留時間減短，氣固接觸效率降低。

圖5給出了流化床內不同水平截面處固相濃度時均值的徑向分布。從圖中可以看出，沿床層徑向，中心區(qū)固含率較壁面附近小，顆粒濃度呈現(xiàn)出中心稀、邊壁濃的基本流動規(guī)律。沿床層軸向，顆粒濃度顯示出上稀下濃的現(xiàn)象，床層下部區(qū)域由于氣泡的存在，顆粒濃度徑向變化更為劇烈，床層中、上部由于氣泡的破碎，固相濃度分布更加均勻。隨著流化床長度的增加，由于氣泡數(shù)量的增多，床層下部固含率波動變大。

3.3 固相速度分布

圖4 t=1000 s時不同尺寸流化床中心豎直截面處固相濃度分布Fig.4 Contour of solid volume fraction at t=1000 s in central vertical plane for fluidized beds with different scale

圖5 不同尺寸流化床的固相濃度徑向分布Fig.5 Radial distribution of solid volume fraction in fluidized beds with different scale

圖6是流化床內不同水平截面處顆粒速度時均值的徑向分布。從圖中可以看出，當流化床床層高度小于（等于）0.15 m 時，長度為0.07 m 和0.15 m 的流化床床層中心區(qū)顆粒向上運動，邊壁區(qū)顆粒向下運動，并且邊壁區(qū)較寬；而對于長度為0.31 m 和0.63 m 的流化床，在較窄的邊壁區(qū)，顆粒向下運動，在中心區(qū)與邊壁區(qū)之間有速度波峰出現(xiàn)，在中心區(qū)部分顆粒被氣泡夾帶而向上運動，在氣泡兩側及下部顆粒向下流動。當流化床高度大于（等于）0.2 m時，四個不同尺寸流化床內所有顆粒均向下流動，這是因為在床層上部氣泡破碎，顆粒只在重力的作用下自由下落，而當接近床層表面（z=0.25 m）時，速度徑向分布則變得較為均勻。

圖6 不同尺寸流化床的固相速度徑向分布Fig.6 Radial distribution of solid velocity in fluidized beds with different scale

3.4 顆粒停留時間分布

實驗中采用脈沖法統(tǒng)計顆粒停留時間分布時，需要將示蹤劑在很短的時間內從流化床入口注入系統(tǒng)。而在CFD-DEM 模擬中，顆粒本身充當示蹤劑，由于入口處要以給定的速度生成顆粒，因此在很短的時間內如果在入口插入大量顆粒將會改變顆粒的進料速率。此外，一次性生成太多顆粒也會對流化床的流體力學行為產生很大影響。因此，在實際模擬中，當系統(tǒng)運行達到穩(wěn)態(tài)后，把在一定時間間隔（150 s）內從入口生成的顆粒標記為示蹤顆粒，然后監(jiān)控這些顆粒的運動與停留時間，將其用于顆粒RTD的分析。

圖7 是不同大小流化床內顆粒停留時間分布，E(t)為停留時間分布密度函數(shù)(s-1)，采用式（32）計算。

式中，mtracer,i為30 s 內從出口流出的示蹤顆粒i的質量;n 為流出示蹤顆粒數(shù)量;mtracer為該時間間隔內進入系統(tǒng)的示蹤顆粒總質量。可以看出，四個不同尺寸的流化床顆粒RTD 出峰都很早，出現(xiàn)長拖尾現(xiàn)象，流動向全混流靠近，反映出典型的鼓泡流化床所具有的RTD特征[23]。在其他條件保持不變的情況下，隨著流化床長度的增加，顆粒RTD 峰值降低，分布變寬。一方面由于床內細顆粒返混嚴重[20,39]，另一方面粗顆粒所受氣體曳力較小，從而導致停留時間更長，曲線長拖尾明顯。這意味著顆粒逆流嚴重而不能及時流出系統(tǒng)，對于含有化學反應的多分散流動系統(tǒng)，可能會造成粗顆粒過度反應而細顆粒轉化不完全，從而降低反應效率。此外，通過與實驗數(shù)據(jù)[圖7（b）]的對比，發(fā)現(xiàn)所提出的Di Felice-Ganser曳力模型很好地預測了0.15 m流化床的顆粒RTD。實驗中第一個數(shù)據(jù)點出現(xiàn)異常，這是因為流動穩(wěn)定后大量示蹤劑的注入引起床層波動所致，而在模擬中所采用的示蹤顆粒進料方式并不會給顆粒整體的流動帶來任何影響。從RTD 曲線的光滑程度來看，流化床越大，曲線噪聲越多，這是因為較大的流化床示蹤顆粒數(shù)量不足，即在流化床尺寸放大的同時，示蹤顆粒的質量（數(shù)量）同樣需要“放大”。在之前的研究中發(fā)現(xiàn)[34]，示蹤顆粒進料時間間隔越長，示蹤顆粒數(shù)量越多（即用于計算停留時間分布密度函數(shù)的顆粒樣本數(shù)越多），RTD 曲線越光滑，但在一定時間段內，示蹤顆粒數(shù)量對曲線的峰值、寬度以及顆粒平均停留時間幾乎沒有影響。因此，如有必要可以通過標記更多的示蹤顆粒來減少較大流化床顆粒RTD曲線的噪聲。

圖7 不同尺寸流化床顆粒停留時間分布Fig.7 Particle RTD of fluidized beds with different scale

3.5 平均停留時間

表4列出了流化床放大時顆粒平均停留時間的模擬結果。本文中平均停留時間tm采用式（33）計算。

表4 不同尺寸流化床顆粒平均停留時間Table 4 Particle MRT of fluidized beds with different scale

圖8 不同粒徑顆粒MRT與流化床長度的關系Fig.8 Relationship between MRT of particles with different sizes and bed length

為進一步發(fā)現(xiàn)放大過程中不同尺寸流化床顆粒MRT 之間的關系，圖8 給出了寬篩分系統(tǒng)每種顆粒以及所有顆粒MRT 隨流化床長度的變化曲線。從圖中可以看出，無論每種顆粒還是所有顆粒，其平均停留時間與流化床長度均表現(xiàn)出線性關系，這是由于固體顆粒進料速率保持不變，顆粒橫向平均速度也基本不發(fā)生改變，系統(tǒng)內每種顆粒以及所有顆粒平均停留時間便隨流化床長度線性增加。可以利用這個結果來預測更大尺寸同類型流化床顆粒的平均停留時間。另外可以看到，流化床越長，不同粒徑顆粒MRT 的差異越大，說明流化床長度的變化對顆粒停留時間也起到了一定的調控作用。

4 結 論

采用耦合多分散、非球形顆粒曳力模型的粗粒化CFD-DEM 方法研究了不同尺寸（長度）的連續(xù)操作流化床流體力學特性和顆粒停留時間及其分布，從中可以得出以下結論。

（1）通過與實驗數(shù)據(jù)的對比，發(fā)現(xiàn)Di Felice-Ganser 曳力模型能夠較為準確地預測多分散、非球形顆粒系統(tǒng)的流動行為和停留時間分布，全床壓降和平均停留時間與實驗結果的相對誤差分別為6.7%和4.6%。

（2）通過顆粒濃度徑向分布的模擬結果可以發(fā)現(xiàn)，流化床越長，顆粒濃度隨徑向位置的波動越劇烈，非均勻結構越明顯。

（3）從不同尺寸流化床顆粒停留時間分布及平均停留時間可以看出，隨著流化床長度的增加，顆粒返混增強，MRT 增大，流化床長度與MRT 呈線性關系。對于同一系統(tǒng)中不同粒徑的顆粒，其MRT 之間的差異隨流化床的放大而增大。

符 號 說 明

D——流化床水力直徑，m

dCGP——粗顆粒直徑，m

Fc——接觸力，N

Fdrag——曳力，N

ICGP,i——轉動慣量，kg·m2

p——氣相壓力，Pa

Rep——顆粒Reynolds數(shù)

tm——顆粒平均停留時間，s

ug——表觀氣速，m/s

umf——最小流化速度，m/s

εmf——最小流化空隙率

μg——氣相黏度，Pa·s

ρg——氣相密度，kg/m3

ρp——顆粒真實密度，kg/m3

τg——氣相應力張量，Pa

下角標

CGP——粗顆粒

g——氣相

p——真實顆粒