基于雙歐拉模型的流化床氣固兩相流場(chǎng)數(shù)值模擬

朱 博， 王 強(qiáng)， 郭慧軍， 張瀟玲， 杜明俊

( 1. 遼寧石油化工大學(xué) 化學(xué)與材料科學(xué)學(xué)院，遼寧 撫順 113001； 2. 中國(guó)石油集團(tuán)工程設(shè)計(jì)有限責(zé)任公司 華北分公司，河北 任丘 062552 )

0 引言

目前,對(duì)流化床床層內(nèi)復(fù)雜動(dòng)力學(xué)特性的研究，雖然取得一定的進(jìn)展，但存在許多問題，如氣固相間曳力的選擇、顆粒相間作用機(jī)制、邊界條件處理、湍流模型選取等，研究大多依賴經(jīng)驗(yàn)技術(shù)和半經(jīng)驗(yàn)理論.隨著計(jì)算機(jī)科學(xué)及數(shù)值分析的不斷發(fā)展，逐漸形成氣固兩相流動(dòng)數(shù)值模擬學(xué)科，預(yù)報(bào)流化床工作的全過程，并為工程流態(tài)化裝置的優(yōu)化和放大設(shè)計(jì)提供一種新的手段.

人們對(duì)氣固流化床內(nèi)兩相流場(chǎng)進(jìn)行研究[1-8]，并取得一些有意義的成果.Alder B J等[4]首次在分子系統(tǒng)相間擴(kuò)散數(shù)值模擬中提出硬球模型.Moon S J等[5]采用數(shù)值方法研究三維振動(dòng)流化床中黏性顆粒的動(dòng)力學(xué)行為,并發(fā)現(xiàn)瞬時(shí)振動(dòng)產(chǎn)生的巨大張力有助于降低黏性顆粒形成的聚團(tuán).何成銓[6-7]等研究細(xì)粉催化劑在流化床中的密度分布規(guī)律，并得出分布圖和關(guān)系式.Gidaspow D[8]提出把顆粒作為連續(xù)相處理的雙流體模型，并且成功預(yù)測(cè)鼓泡床流體力學(xué)行為.楊太陽(yáng)等[9]在雙流體模型基礎(chǔ)上發(fā)展針對(duì)氣泡的模擬程序，成功模擬鼓泡床中氣泡頻率、速度、長(zhǎng)度特征及顆粒循環(huán)狀態(tài).

圖1 流化床物理模型(單位:m)

對(duì)于氣固兩相流動(dòng)，相間曳力是表征兩相之間相互作用和動(dòng)量交換的重要參數(shù)，決定數(shù)值計(jì)算總體的準(zhǔn)確性[10].對(duì)于低流化風(fēng)速下，采用Syamlal-O’Brien曳力模型計(jì)算床內(nèi)流化過程誤差較大[11-12].結(jié)合實(shí)驗(yàn)數(shù)據(jù)，筆者對(duì)Syamlal-O’Brien曳力模型進(jìn)行修正，利用FLUENT軟件自帶的雙歐拉模型對(duì)均勻床內(nèi)低速氣固兩相流場(chǎng)進(jìn)行數(shù)值計(jì)算，為進(jìn)一步改善和推動(dòng)氣固兩相流動(dòng)數(shù)值模擬的發(fā)展提供依據(jù).

1 模型修正

實(shí)驗(yàn)和數(shù)值模擬采用高1.00 m，直徑為0.15 m帶分布板的鼓泡床，入口空氣流速為0.25 m/s.頂部壓力為標(biāo)準(zhǔn)大氣壓，固相區(qū)充滿粒狀固體，體積分?jǐn)?shù)為55%，簡(jiǎn)化的流化床物理模型見圖1.利用Gambit進(jìn)行建模和網(wǎng)格劃分.采用結(jié)構(gòu)化矩形網(wǎng)格進(jìn)行單元?jiǎng)澐?，并通過床層壓降確定模型最適宜的網(wǎng)格數(shù)(30 480個(gè)).

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

應(yīng)用較廣泛的氣固兩相流模型包括歐拉-歐拉擬流體和歐拉-拉格朗日離散相模型，基于計(jì)算機(jī)資源的限制，大多采用歐拉-歐拉擬流體模型進(jìn)行數(shù)值計(jì)算[13-14].在實(shí)際流態(tài)化過程中，薩夫曼升力對(duì)于大粒徑顆粒影響較大，對(duì)于粒徑遠(yuǎn)小于粒子間距的鼓泡床，可忽略升力對(duì)主相流場(chǎng)的影響[14].對(duì)于虛擬質(zhì)量力，當(dāng)?shù)诙嗝芏冗h(yuǎn)小于主相密度時(shí)，虛擬質(zhì)量影響較大，不能忽略，鼓泡床中第二相顆粒的密度遠(yuǎn)大于主相氣體的密度，因此該項(xiàng)可以忽略[14].采用雙歐拉模型對(duì)低流化風(fēng)速下氣固兩相流場(chǎng)進(jìn)行數(shù)值計(jì)算，可以忽略升力、虛擬質(zhì)量力及壁面滑移邊界對(duì)流態(tài)化過程的影響.

質(zhì)量守恒方程為

(1)

動(dòng)量守恒方程為

(2)

(3)

式中：αv，αs分別為氣相和固相的體積分?jǐn)?shù)；uv,us分別為氣相和固相的速度矢量；ρv，ρs分別為氣相和固相的密度；p為兩相共享壓力；τ為壓力應(yīng)變張量；Rvs，Rsv分別為氣相和固相間曳力；分別為氣相和固相間質(zhì)量傳遞；usv，uvs分別為氣相和固相間速度；g為重力加速度.

1.2 修正Syamlal-O’Brien阻力模型

最小流化速度是指流化床層開始流化時(shí)的氣體空床線速.原Syamlal-O’Brien阻力模型計(jì)算的最小流化速度為0.21 m/s，不能滿足低速流化床相間阻力的計(jì)算，結(jié)合實(shí)驗(yàn)數(shù)據(jù)，重新對(duì)相對(duì)雷諾數(shù)和曳力系數(shù)進(jìn)行界定，對(duì)原有的Syamlal-O’Brien阻力模型進(jìn)行修正，使其滿足低速流化過程.

曳力f計(jì)算模型為

(4)

式中：CD為曳力函數(shù)；Res為相對(duì)雷諾數(shù)；vr,s為固相沉降速率.

曳力函數(shù)CD采用Dalla Valle給出的形式：

(5)

這個(gè)模型是基于流化床或沉淀床顆粒的末端速度的測(cè)量，并使用體積分?jǐn)?shù)和相對(duì)雷諾數(shù)的函數(shù)關(guān)系式：

(6)

式中：ρl為液相黏度；ds為固相顆粒直徑；vs為固相速度；vl為液相速度；μl為液相黏度.

液—固交換系數(shù)Ksl為

(7)

式中：αs,αl分別為固相和液相分?jǐn)?shù).

固相自由沉降速度系數(shù)為

(8)

1.3 邊界條件

計(jì)算時(shí)，固相顆粒密度為2 600 kg/m3,剪切黏度為1.789×10-5Pa·s,粒徑為0.2 mm,氣相密度為1.2 kg/m3，剪切黏度為1.89×10-5Pa·s.

入口邊界條件為

出口邊界條件為

2 數(shù)值模擬

實(shí)驗(yàn)裝置由流化床和輔助設(shè)備組成，靜床高度為1.00 m,直徑為0.15 m,固體顆粒采用Geldart B類顆粒，空氣為流化氣體，實(shí)驗(yàn)在常溫常壓下進(jìn)行，實(shí)驗(yàn)裝置及流程見圖2.通過實(shí)驗(yàn)測(cè)得最小流化速度為0.80 m/s.

數(shù)值計(jì)算采用FLUENT 6.3軟件，并對(duì)其進(jìn)行二次開發(fā)，利用C語(yǔ)言編寫UDF程序?qū)崿F(xiàn)修正后Syamlal-O’Brien模型的計(jì)算過程，采用一階迎風(fēng)格式對(duì)控制方程進(jìn)行離散，利用Phass Coupled SIMPLE算法進(jìn)行求解，迭代步長(zhǎng)為1 ms,修正壓力和速度的松弛因子，使其快速收斂.

穩(wěn)定進(jìn)氣10 s后流化床軸向壓力變化的實(shí)驗(yàn)結(jié)果與計(jì)算結(jié)果見圖3.由圖3可知，采用修正Syamlal-O’Brien曳力模型計(jì)算結(jié)果小于實(shí)驗(yàn)結(jié)果，但軸向壓降變化趨勢(shì)相同，這主要是受固體顆粒球形度的影響.對(duì)比計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果，最大相對(duì)誤差為4.4%，滿足工程實(shí)際要求，修正模型能夠描述低速氣固兩相流動(dòng)相互作用物理過程.

圖2 流化床實(shí)驗(yàn)裝置及流程

圖3 10 s后流化床軸向壓力變化曲線

Syamlal-O’Brien阻力模型修正前后不同時(shí)刻流化床固相顆粒體積濃度云圖見圖4.由圖4可知，修正前模型在0.8 s后開始出現(xiàn)氣泡，且隨著時(shí)間延長(zhǎng)，氣泡不斷上升并長(zhǎng)大；修正后模型在1.4 s后逐漸產(chǎn)生氣泡，且氣泡增長(zhǎng)和上升速度較慢.對(duì)于固相顆粒，修正前模型在1.5 s前均勻膨脹，之后在床底不均勻氣流的脈動(dòng)下，開始出現(xiàn)不均勻膨脹，且密相床層高度隨之降低；修正后模型在2.0 s前上部顆粒均勻膨脹.流態(tài)化初期，密相床層軸向上存在速度差，導(dǎo)致高度方向上流態(tài)化氣速大小不一，較低的氣速使流化床床層發(fā)生均勻膨脹，較高的氣速導(dǎo)致氣泡生成.由于計(jì)算所用顆粒直徑為0.2 mm，屬于流態(tài)化原理中Geldart B類顆粒，該類顆粒在流態(tài)化氣速較小時(shí)流化床床層先發(fā)生均勻膨脹而無氣泡生成，隨著流態(tài)化氣速增大至鼓泡速度時(shí)才出現(xiàn)氣泡.修正后模型計(jì)算結(jié)果與流態(tài)化原理[15]導(dǎo)出的結(jié)論一致.

穩(wěn)定進(jìn)氣2.0 s后不同介質(zhì)流化床軸向壓力變化曲線見圖5.由圖5可知，靠近床層底部一定區(qū)域，修正后的曳力模型計(jì)算的固相和氣相床層壓力大于原模型的，且隨床層高度的增加壓力逐漸降低，主要是修正曳力模型計(jì)算的氣泡聚并與破裂頻率降低引起的.流化床床內(nèi)壓力場(chǎng)的分布受氣泡狀態(tài)及進(jìn)氣氣流影響，采用原曳力模型計(jì)算的固相介質(zhì)和氣相介質(zhì)的軸向壓力波動(dòng)較大，主要原因是修正前模型計(jì)算的最小流化速度偏大，短時(shí)間內(nèi)出現(xiàn)大量氣泡，隨著氣泡的不斷上升、融合和破碎，使固相區(qū)出現(xiàn)不均勻膨脹，從而導(dǎo)致床層壓力開始大幅度不規(guī)則波動(dòng)，此時(shí)床層壓降受脈動(dòng)氣流頻率和氣泡行為共同控制.修正后模型計(jì)算壓力波動(dòng)的范圍較小，主要原因是初始流化速度小，床層壓力主要受脈動(dòng)氣流波動(dòng)控制，使床內(nèi)固相顆粒均勻膨脹.

圖4 不同時(shí)刻流化床固相顆粒體積濃度云圖

圖5 2.0 s后不同介質(zhì)流化床軸向壓力變化曲線

3 結(jié)論

(1)結(jié)合實(shí)驗(yàn)數(shù)據(jù)，對(duì)體積分?jǐn)?shù)和相對(duì)雷諾數(shù)進(jìn)行重新界定，修正Syamlal-O’Brien曳力模型，使其最小流化速度為0.80 m/s,滿足低流化風(fēng)速條件.

(2)結(jié)合修正后曳力模型，采用雙歐拉模型對(duì)均勻床內(nèi)低速氣固兩相流場(chǎng)進(jìn)行數(shù)值計(jì)算，計(jì)算結(jié)果比實(shí)驗(yàn)數(shù)據(jù)略小，最大相對(duì)誤差為4.4%，滿足工程實(shí)際要求.

(3)由于Syamlal-O’Brien阻力模型是根據(jù)床層壓降、表觀氣速、床層膨脹比等實(shí)驗(yàn)數(shù)據(jù)得到的半經(jīng)驗(yàn)公式，是在單個(gè)顆粒曳力模型的基礎(chǔ)上引入顆粒體積分?jǐn)?shù)描述對(duì)周圍顆粒的影響，建議結(jié)合實(shí)驗(yàn)，根據(jù)數(shù)學(xué)理論和數(shù)值模擬等繼續(xù)完善一個(gè)相對(duì)完整的阻力模型.