氣-液-固三相流混合建模與求解方法*

范興華 譚大鵬 李霖 殷梓超 王彤

(浙江工業(yè)大學(xué)機(jī)械工程學(xué)院, 特種裝備制造與先進(jìn)加工技術(shù)教育部/浙江省重點(diǎn)實(shí)驗(yàn)室, 杭州 310014)

1 引 言

氣-液-固三相流混合是高端化工、鋰電生產(chǎn)的關(guān)鍵工藝環(huán)節(jié), 混合執(zhí)行構(gòu)件和流道物理空間需要提供較高的傳質(zhì)和高湍流能力, 且伴隨強(qiáng)剪切過程.上述要求使得氣-液-固三相流混合過程非常復(fù)雜, 且難于觀察整體流場及關(guān)鍵區(qū)域的顆粒分布[1-3].混合物理空間幾何尺度相對(duì)于顆粒要高多個(gè)數(shù)量級(jí), 其內(nèi)部三維循環(huán)流動(dòng)和湍流多相流的復(fù)雜性, 給混合執(zhí)行機(jī)構(gòu)優(yōu)化設(shè)計(jì)、混合過程邊界條件調(diào)控提出了重要技術(shù)挑戰(zhàn)[4,5].

當(dāng)前計(jì)算流體力學(xué)(computational fluid dynamics, CFD)方法廣泛應(yīng)用于液-固混合過程的模擬計(jì)算, 一般基于歐拉-歐拉模型, 將顆粒固相視為連續(xù)相, 來描述相間相互滲透過程, 該模型占用計(jì)算資源比較少, 但模擬精度較低, 且無法獲得顆粒運(yùn)動(dòng)狀態(tài)[6-10].基于歐拉-拉格朗日模型的離散單元法(discrete element method, DEM)可以獲得顆粒的運(yùn)動(dòng)和相互作用, 可以與不同的流體動(dòng)力學(xué)計(jì)算方法相結(jié)合, 來模擬流體-顆粒流[11,12], 如格子-玻爾茲曼方法(lattice Boltzmann method, LBM).該方法在離散的晶格網(wǎng)格上使用代表流體相的虛擬顆粒, 并通過求解離散的Boltzmann方程模擬流體的流動(dòng)[13].相關(guān)學(xué)者已經(jīng)對(duì)三維LBM-DEM耦合進(jìn)行了嘗試[14], 但是由于要求流體尺寸要比固體顆粒尺寸小得多, 對(duì)計(jì)算能力要求非常高, 多數(shù)針對(duì)三維問題的LBM-DEM耦合解法仍在開發(fā)中.將CFD與離散單元法(CFD-DEM)耦合使用,可以預(yù)測顆粒尺度的變化, 顆粒-顆粒和顆粒-壁之間的相互作用都通過牛頓運(yùn)動(dòng)方程求解, 而顆粒-流體之間的相互作用則通過源相交換來實(shí)現(xiàn).Bastien等[15]使用CFD-DEM模型, 以非常好的可靠性再現(xiàn)了固-液混合系統(tǒng)中發(fā)生的各種現(xiàn)象,研究證明, 在非慣性參考系下進(jìn)行CFD-DEM模擬是可行的, 這為CFD-DEM的應(yīng)用開辟了廣闊的前景.Shao等[16]利用CFD-DEM耦合的模型研究了三維混合過程中的固體懸浮行為, 與實(shí)驗(yàn)測量和基于歐拉-歐拉方案的CFD仿真相比, CFDDEM仿真可以提供更多流場信息.Blais等[17-19]開發(fā)了一種半解析CFD-DEM模型, 并利用該模型進(jìn)行固-液混合操作的設(shè)計(jì)以及優(yōu)化, 提高了顆粒懸浮比例, 改善了流型和顆粒分布, 但未考慮自由液體表面的分布及其穩(wěn)定性.

對(duì)于自由液體表面流, 常用流體體積(volume of fluid, VOF)模型進(jìn)行模擬, Xu等[20]研究了表面渦流對(duì)混合物理空間中顆粒分散的影響, 研究發(fā)現(xiàn)表面渦流的產(chǎn)生降低了顆粒分布的均勻性, 混合擋流物理構(gòu)件可使得表面渦流得到有效控制, 提高了顆粒分布的均勻性.Sun和Sakai[21]開發(fā)了一種基于歐拉-拉格朗日方法的模型, 可以模擬復(fù)雜的三相流行為, 包括自由表面的變形和顆粒引起的液體位移.Wu等[22]利用虛擬雙重網(wǎng)格孔隙度模型改進(jìn)了DEM-VOF模型, 新模型克服了計(jì)算過程中的不穩(wěn)定性, 可用于固-液混合過程的流場模擬計(jì)算.Kang等[23]利用DEM-VOF模型結(jié)合雷諾壓力模型(RSM), 對(duì)具有自由表面的無擋流構(gòu)件的混合物理空間中的顆粒懸浮動(dòng)力學(xué)進(jìn)行了模擬,得到了混合流道的幾何形狀、葉輪轉(zhuǎn)速、顆粒密度和直徑等對(duì)自由表面渦流、流型和顆粒懸浮動(dòng)力學(xué)的影響規(guī)律.

當(dāng)前關(guān)于多相流的研究主要以兩相為主, 考慮三相耦合過程的多集中于自由液面, 卻鮮有考慮內(nèi)部充氣對(duì)液-固兩相的影響, 尤其是視固相為離散顆粒相的情況.因此, 建立流體和顆粒的動(dòng)力學(xué)模型, 通過求解動(dòng)量方程, 實(shí)現(xiàn)流體與顆粒的雙向耦合, 進(jìn)行氣-液-固三相流混合的研究, 揭示三相耦合在復(fù)雜混合過程中的作用規(guī)律非常必要.

針對(duì)上述目標(biāo), 本文首先建立了氣-液-固三相流動(dòng)力學(xué)模型, 分別包括VOF模型、DEM模型以及兩者的耦合模型; 然后進(jìn)行了建模與網(wǎng)格劃分以及邊界條件和參數(shù)設(shè)置, 并進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證, 進(jìn)而選取最終使用的網(wǎng)格數(shù)量并開始不同案例的數(shù)值計(jì)算; 自主開發(fā)用戶自定義函數(shù)(user defined function, UDF)通信接口, 得到流體與顆粒間的相互作用力, 提出了一種多孔相間耦合解法來描述顆粒運(yùn)動(dòng)軌跡; 最后通過計(jì)算結(jié)果討論了充氣對(duì)自由液面、流體速度以及顆粒懸浮的影響, 揭示了不同充氣條件下混合物理空間內(nèi)多相流的演變規(guī)律并得出了相應(yīng)結(jié)論.

2 三相流動(dòng)力學(xué)模型

對(duì)于VOF-DEM模型耦合的多相流, 通過對(duì)體積平均的Navier-Stokes方程進(jìn)行求解, 進(jìn)而對(duì)連續(xù)相的流體進(jìn)行描述, 針對(duì)VOF模型還存在自由液面的問題, 氣-液兩相存在明顯交界面, 可通過2.1節(jié)模型求解氣-液兩相問題.同時(shí)使用離散單元法對(duì)固體相顆粒進(jìn)行建模, 可通過2.2節(jié)模型求解.這兩個(gè)模型以一定的時(shí)間步間隔進(jìn)行雙向耦合交換數(shù)據(jù), 一般CFD的時(shí)間步長明顯大于DEM時(shí)間步長, 以便正確獲得接觸作用, 顆粒通常不會(huì)在單個(gè)DEM時(shí)間步長中移動(dòng)很遠(yuǎn)的距離.因此, 不需要兩者1∶1的時(shí)間步長比, 對(duì)于DEM,CFD的時(shí)間步長比一般從1∶10到1∶100不等, 兩者耦合時(shí)選擇合適的時(shí)間步長比, 不僅可以節(jié)省計(jì)算時(shí)間, 還可以避免計(jì)算的發(fā)散.

2.1 VOF模型

混合過程中存在復(fù)雜的氣-液兩相耦合現(xiàn)象,所以應(yīng)該用多相流模型來描述.VOF模型是基于歐拉網(wǎng)格下的表面追蹤模型, 通過求解單相或多相的體積分?jǐn)?shù)來追蹤和捕捉互不相融流體間的交界面.通過計(jì)算各相的控制方程, 能夠準(zhǔn)確地模擬混合過程中多相組分的動(dòng)態(tài)演化和瞬態(tài)特征的捕捉,其中流體的連續(xù)性方程和動(dòng)量方程表示如下[21]:

式中,ρf是流體密度,εf是空隙率,u是流體速度,p是壓力,μ是流體動(dòng)力黏度,Fpf是顆粒流體間的相互作用項(xiàng)的反作用力,Fst是自由液面附近的表面張力.為了提高模擬精度, 使計(jì)算更加接近實(shí)際情況, 本文采用連續(xù)表面張力(continuum surface force, CSF)模型處理表面張力, 其表達(dá)式如下:

式中σ是流體的表面張力系數(shù); 而κ是氣-液相界面的曲率, 其表達(dá)式為

其中n=?·α2是次相體積分?jǐn)?shù)的法向量.

對(duì)于VOF模型氣-液兩相的交界面可通過界面?zhèn)鬏敺匠糖蠼?

其中α是液體的體積分?jǐn)?shù), 若α=1 , 代表該計(jì)算單元全是液相, 不含氣相; 若 0＜α＜1 , 則說明該計(jì)算單元內(nèi)同時(shí)包含氣、液兩相; 若α=0 , 代表該相全部是氣相.因此, 基于VOF模型可以求解氣-液交界面的相互作用過程, 尤其是在強(qiáng)剪切作用下的氣相剝離過程.

此外, 多相流環(huán)境中尤其是強(qiáng)剪切區(qū)域處于湍流狀態(tài)下, 為了能夠得到精確的模擬結(jié)果, 流體控制方程選擇標(biāo)準(zhǔn)的湍動(dòng)能-耗散(k-ε)模型, 該模型在劇烈變化的流場中有較好的計(jì)算性能, 其控制方程如下[24]:

式中,k是湍動(dòng)能;ε是湍動(dòng)能耗散率;Gk是由于平均速度梯度引起的湍動(dòng)能產(chǎn)生項(xiàng);Gb是由于浮力影響引起的湍動(dòng)能產(chǎn)生項(xiàng);YM為可壓縮湍流脈動(dòng)膨脹對(duì)總的耗散率的影響;C1ε,C2ε,C3ε為經(jīng)驗(yàn)常數(shù), 張量下標(biāo)i,j表示x,y,z三個(gè)方向分量;σk,σε為湍動(dòng)能和湍動(dòng)耗散率對(duì)應(yīng)的普朗特?cái)?shù).

2.2 DEM模型

在考慮流體相時(shí), 已經(jīng)引入了顆粒與流體間的相互作用, 因此在該部分將分析顆粒模型.DEM是一種可以用于計(jì)算非連續(xù)顆粒的運(yùn)動(dòng)規(guī)律, 并且可以分析離散顆粒接觸力以及運(yùn)動(dòng)的分析模型.此模型中, 將顆粒相視為離散相, 相比于其他方法更加貼近實(shí)際, 更能還原顆粒的真實(shí)運(yùn)動(dòng)情況, 顆粒的運(yùn)動(dòng)是基于牛頓第二定律的計(jì)算得出的, 通過計(jì)算可以得到顆粒的平移、旋轉(zhuǎn)的速度和位置隨時(shí)間的變化關(guān)系, 顆粒的平移運(yùn)動(dòng)取決于作用在其上的力的總和, 而旋轉(zhuǎn)運(yùn)動(dòng)則由接觸的轉(zhuǎn)矩控制, 其控制方程可表示為[16,25]

式中,mi是顆粒i的質(zhì)量;xi是顆粒的位移;Fc,ij是顆粒i和j之間的接觸力;Fpf,i是顆粒i與流體間的相互作用力;g是重力加速度;Ii是顆粒i的慣性力矩;θp,i是顆粒i的角位移;Tt,ij和Tr,ij是作用在顆粒上的切向和滾動(dòng)摩擦力矩, (8)式和(9)式中包含的一些力和扭矩的詳細(xì)計(jì)算見如下公式.接觸力:

法向接觸力:

切向接觸力:

切向摩擦力矩:

滾動(dòng)摩擦力矩:

顆粒-流體作用力:

這里Y*表示接觸顆粒間等效楊氏模量;R*表示接觸顆粒間等效半徑;δcn,ij,δct,ij分別表示接觸顆粒間法向、切向重疊大小;Scn,ij,Sct,ij分別表示接觸顆粒間法向、切向接觸剛度;m*表示接觸顆粒間等效質(zhì)量;vcn,vct分別表示接觸顆粒間的法向、切向相對(duì)速度;Lij表示接觸顆粒間中心距離;nij表示接觸顆粒間單位矢量;μr表示顆粒的滾動(dòng)摩擦系數(shù);ωij表示顆粒接觸平面上的角速度矢量; ΔV是顆粒i所在的計(jì)算網(wǎng)格的體積;np表示顆粒數(shù)量;fpf,i表示顆粒i受到的所有流體、固體對(duì)顆粒作用之和, 包括曳力fd,i[26,27]、壓力梯度力f?p,i、黏性力f?τ,i、表面張力fst,i、虛擬質(zhì)量力fvm,i、Basset力fB,i、Saffman升力fSaff,i[28]和Magnus升力fMag,i[23].由于動(dòng)量方程表達(dá)式已經(jīng)直接包含了壓力、黏性力和表面張力, 所以在Fpf中相應(yīng)地把這些力減去了.

2.3 耦合模型

為了獲得流體與顆粒間準(zhǔn)確的相互作用力, 提出了一種相間耦合解法—多孔模型來描述精確的顆粒運(yùn)動(dòng)軌跡, 其計(jì)算表達(dá)式如下:

式中,εps,i代表流體單元內(nèi)第i個(gè)多孔球單位體積內(nèi)的顆粒體積.該模型克服了傳統(tǒng)方法中當(dāng)顆粒粒徑接近網(wǎng)格尺寸時(shí)引起的計(jì)算不穩(wěn)定性問題, 改善了顆粒-流體的相互作用, 并且在計(jì)算過程中, 把顆粒的體積考慮在內(nèi), 因此, 用該方法求解出的流場也更加精確.

將2.1節(jié)和2.2節(jié)兩個(gè)模型以及多孔模型的控制方程寫入接口程序, 進(jìn)行編譯, 最終通過動(dòng)量方程中的動(dòng)量源交換項(xiàng)實(shí)現(xiàn)雙向耦合, 這樣VOFDEM耦合通過編好的用戶自定義函數(shù)(UDF)進(jìn)行數(shù)據(jù)通信, 實(shí)現(xiàn)歐拉雙流體相和拉格朗日顆粒相的雙向耦合.

整體模型的計(jì)算流程如上圖1所示.首先對(duì)流場和顆粒場進(jìn)行初始化, 該過程通過CFD計(jì)算接口文件實(shí)現(xiàn); 然后, 開始計(jì)算, 通過2.1節(jié)(1)式和(2)式迭代求解得到流場的速度和力等信息, 求解(5)式得到關(guān)于自由液面的變化信息, 通過(6)式和(7)式計(jì)算流體的湍動(dòng)能-耗散; DEM通過利用2.2節(jié)的(8)式和(9)式迭代計(jì)算獲得顆粒的速度和位置等信息, 并進(jìn)行更新; 接著通過判斷收斂與否, 進(jìn)行選擇, 如果不收斂通過求解2.3節(jié)(17)式得到流體單元的空隙率, 繼續(xù)前述流場的計(jì)算, 如此循環(huán)實(shí)現(xiàn)雙向耦合, 互相交換數(shù)據(jù), 直到收斂停止計(jì)算, 完成模擬.

圖1 VOF-DEM耦合計(jì)算流程圖Fig.1.VOF-DEM coupling calculation flowchart.

3 數(shù)值計(jì)算

3.1 物理模型與網(wǎng)格劃分

研究所選取的物理混合空間為帶擋流板的半圓底容器, 混合執(zhí)行構(gòu)件為葉輪, 結(jié)構(gòu)如圖2所示,具體物理參數(shù)如下: 直徑T= 200 mm, 高度H=3T/2, 液位高度hl=T, 葉輪直徑D=T/2, 槳葉長度L= 45 mm, 寬度W=T/10, 厚度t= 2 mm,傾斜角為45°, 安裝高度C= 93 mm, 攪拌軸直徑d1= 14 mm, 底部進(jìn)氣口直徑d2= 14 mm, 高度ha= 4 mm, 擋流板高度hb= 11T/10, 寬度為T/10.

圖2 混合空間結(jié)構(gòu)示意圖Fig.2.Diagram of mixed space structure.

首先建立混合物理空間三維模型, 對(duì)流體域進(jìn)行網(wǎng)格劃分, 最終生成的網(wǎng)格如圖3所示.為了方便計(jì)算, 流體域被劃分為包含混合葉輪的轉(zhuǎn)子區(qū)域和除此以外的靜子區(qū)域兩部分.然后對(duì)流體域進(jìn)行離散化處理, 由于轉(zhuǎn)子區(qū)域的變化梯度較大, 尤其是混合葉輪這種小尺寸, 強(qiáng)剪切區(qū)域變化梯度最大, 劃分時(shí)要特別注意, 所以要進(jìn)行局部網(wǎng)格的加密處理, 劃分后的靜子區(qū)域網(wǎng)格尺寸為7 mm, 轉(zhuǎn)子區(qū)域網(wǎng)格尺寸為5 mm, 葉面網(wǎng)格尺寸進(jìn)一步細(xì)化為3 mm, 進(jìn)氣口也進(jìn)行適當(dāng)?shù)募用? 網(wǎng)格尺寸為3 mm, 劃分網(wǎng)格后, 網(wǎng)格的正交質(zhì)量保持在0.5以上, 最終用于計(jì)算的網(wǎng)格數(shù)目為31萬.

圖3 網(wǎng)格劃分 (a) 靜子區(qū)域網(wǎng)格; (b) 轉(zhuǎn)子區(qū)域網(wǎng)格;(c) 葉輪網(wǎng)格Fig.3.Grids division: (a) Grids of stator region; (b) grids of rotor region; (c) grids of impeller.

3.2 邊界條件及參數(shù)選擇

混合過程模擬采用了瞬態(tài)的VOF模型計(jì)算,選擇顯式的時(shí)間離散格式, 湍流模型使用標(biāo)準(zhǔn)的湍動(dòng)能-耗散(k-ε)模型, 該模型具有較高的物理可靠性, 可為混合湍流過程提供精確解, 近壁區(qū)域選擇標(biāo)準(zhǔn)壁面函數(shù).混合容器頂部設(shè)置為壓力出口邊界條件, 容器壁面為無滑移壁面邊界條件, 底部為充氣管道, 并采用速度入口邊界條件.數(shù)值求解方法使用coupled方案, 該解法耦合了動(dòng)量、能量及組分方程, 能比較快地得到收斂解, 動(dòng)量離散格式、湍動(dòng)能和湍動(dòng)能耗散率離散格式均采用二階迎風(fēng)以獲得精確的解[29-31], 體積分?jǐn)?shù)離散格式使用分段線性界面重構(gòu)(piecewise linear interface construction, PLIC)算法, 這種方法是精度最高的一種[24,29], 監(jiān)視器收斂殘差均為10—6.

葉輪旋轉(zhuǎn)模型常用到滑移網(wǎng)格(sliding-grid,SG)方法和多重參考系(multiple reference frames,MRF)方法.其中, SG方法常用于瞬態(tài)模擬, 而MRF方法通常用于穩(wěn)態(tài)模擬, 文獻(xiàn)[32]采用兩種方法對(duì)比得到的最終結(jié)果非常相似.MRF方法也能夠用于瞬態(tài)仿真, 此種情況是以偽穩(wěn)態(tài)方式進(jìn)行計(jì)算, 與更準(zhǔn)確、但更耗時(shí)的SG方法相比, 可節(jié)省大量計(jì)算機(jī)資源, 精度能滿足多數(shù)場景的需要[33-35].因此本次模擬采用瞬態(tài)MRF方法進(jìn)行.對(duì)于MRF方法, 需要將流體區(qū)域劃分為內(nèi)部動(dòng)區(qū)域和外部靜止區(qū)域兩部分, 兩部分通過交界面(interface)進(jìn)行數(shù)據(jù)傳遞.對(duì)于本研究所涉及的其他物理參數(shù)設(shè)置, 見表1.

表1 流體和顆粒特性設(shè)置Table 1.Characteristics settings of fluid and particle.

3.3 網(wǎng)格無關(guān)性分析

網(wǎng)格的大小會(huì)直接影響數(shù)值模擬的結(jié)果, 一般情況下, 網(wǎng)格越小, 計(jì)算結(jié)果也越精確, 但是, 隨之帶來的是網(wǎng)格數(shù)量也越來越多, 需要更多的計(jì)算時(shí)間才能收斂, 給計(jì)算帶來了很大的困難.因此, 有必要進(jìn)行網(wǎng)格獨(dú)立性研究, 以確保計(jì)算誤差在可接受的范圍內(nèi), 得到準(zhǔn)確的模擬結(jié)果.本次以傾角45°的槳式葉輪進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證, 模擬轉(zhuǎn)速為400 r/min時(shí)的混合流場, 使用了四種網(wǎng)格尺寸,總網(wǎng)格數(shù)分別為219000, 311000, 425000, 661000,考察網(wǎng)格數(shù)量對(duì)模擬結(jié)果的影響, 對(duì)5 s時(shí)槽內(nèi)Z= 150 mm,Y= 0 mm,X從—100到100 mm的軸向速度場進(jìn)行對(duì)比, 結(jié)果如圖4所示.

圖4 四種網(wǎng)格尺寸在t = 5 s時(shí)的軸向速度分布Fig.4.Axial velocity distribution of four grid sizes at t = 5 s.

可以發(fā)現(xiàn)流場在不同位置的軸向速度具有相似的趨勢, 但網(wǎng)格數(shù)量為219000時(shí)整體具有較大的誤差, 而其他三種網(wǎng)格計(jì)算誤差在5%以內(nèi), 滿足網(wǎng)格獨(dú)立性要求.基于上述結(jié)果, 后面計(jì)算模型所采用的網(wǎng)格數(shù)量均為311000, 這樣既減少了計(jì)算量, 又可以得到比較準(zhǔn)確的結(jié)果.

4 結(jié)果與討論

如前所述, 氣-液-固三相混合物理空間內(nèi)部是一個(gè)復(fù)雜的湍流環(huán)境, 擋流構(gòu)件的存在增加了攪拌速度下湍流場的無序性和非線性, 為了獲得其中的流體-固體多相耦合和相間傳質(zhì)特性, 將深入研究對(duì)比不同充氣條件下對(duì)混合空間內(nèi)自由表面、流體流動(dòng)和固體顆粒懸浮的影響.在數(shù)值算例中, 顆粒直徑均為1 mm, 顆粒數(shù)目均為10000, 且顆粒隨機(jī)分布在混合容器的底部區(qū)域內(nèi), 在初始條件下僅受重力作用.

4.1 對(duì)自由液面的影響

基于數(shù)值模擬的結(jié)果, 首先研究充氣狀態(tài)對(duì)自由液面的影響, 圖5為t= 5 s時(shí)四種充氣條件下混合空間內(nèi)的實(shí)際自由液面圖, 藍(lán)色為自由液面的演化形態(tài), 分別對(duì)應(yīng)充氣速度為0, 0.01, 0.05和1 m/s.為了能夠清晰地看出自由液面的變化, 隱藏了底部液體和顆粒, 只保留了自由液面.

通過圖5(a)—(d)可以看出, 在圖5(a)未充氣的混合空間中液面存在小幅的波動(dòng), 因?yàn)榈撞繑嚢鑼?duì)流場的擾動(dòng)能量向上傳輸至自由液面時(shí), 能量的衰減不足以對(duì)液面造成較大的擾動(dòng), 而這種輕微的波動(dòng)主要是顆粒和液體的相互作用造成的; 對(duì)比之下, 充氣速度v= 0.01和v= 0.05 m/s時(shí), 如圖5(b)和圖5(c)時(shí)兩者液面波動(dòng)都很小, 顯然傳輸?shù)耐膭?dòng)能依然達(dá)不到自由液面的擾動(dòng)閾值; 而當(dāng)充氣速度v= 1 m/s時(shí), 其自由液面如圖5(d)所示, 可以看出液面波動(dòng)非常大, 尤其是擋流板之間的區(qū)域, 液面上升明顯, 有漩渦的產(chǎn)生, 這主要是因?yàn)槌祟w粒對(duì)液面的影響之外, 底部充氣速度較大, 增加了湍流場的流體上沖動(dòng)能, 對(duì)三相流系統(tǒng)造成了較大的擾動(dòng), 且氣體上浮溢出造成了液面的不穩(wěn)定, 進(jìn)而有較大的振蕩.

圖5 t = 5 s時(shí)四種充氣條件下混合空間內(nèi)的自由液面 (a) v = 0 m/s; (b) v = 0.01 m/s; (c) v = 0.05 m/s;(d) v = 1 m/sFig.5.Free surface under four aeration conditions at t =5 s: (a) v = 0 m/s; (b) v = 0.01 m/s; (c) v = 0.05 m/s;(d) v = 1 m/s.

針對(duì)上述充氣擾動(dòng)液面的現(xiàn)象, 考慮與攪拌下流場的流動(dòng)模式密切相關(guān), 圖6給出了不同充氣條件下的切向速度矢量圖.從圖6可以看出, 四種條件下?lián)趿鳂?gòu)件附近的切向速度都是該截面內(nèi)最小的, 是因?yàn)閾趿鳂?gòu)件將流體的切向速度轉(zhuǎn)換為了軸向速度.而當(dāng)高速旋轉(zhuǎn)的流體與擋流構(gòu)件接觸時(shí),快速接觸過程必然引起局部湍流渦團(tuán)的耗散, 增加了流動(dòng)模式的隨機(jī)性, 故擋流構(gòu)件附近的液面振蕩相對(duì)明顯.同時(shí), 通過對(duì)比可以發(fā)現(xiàn),v= 1 m/s時(shí)的流體切向速度分布極不均勻并且最大, 較強(qiáng)的充氣強(qiáng)度對(duì)底部分布的顆粒相起到擾動(dòng)作用, 在攪拌速度和底吹作用下, 整個(gè)流場處于非線性湍流狀態(tài), 這也是引起自由液面漩渦和波動(dòng)的主要原因; 而其他三種狀態(tài)下的切向速度則相對(duì)較小, 其中v= 0.05 m/s時(shí)的切向速度最均勻, 相對(duì)比較穩(wěn)定, 其對(duì)應(yīng)的自由液面也是最穩(wěn)定的, 這也與文獻(xiàn)[23]的結(jié)論相吻合.

圖6 t = 5 s時(shí), 四種充氣條件下混合空間內(nèi)z = 0.15 m高度截面的切向速度矢量圖 (a) v = 0 m/s; (b) v = 0.01 m/s;(c) v = 0.05 m/s; (d) v = 1 m/sFig.6.Tangential velocity vector in height z = 0.15 m under four aeration conditions at t = 5 s: (a) v = 0 m/s; (b) v =0.01 m/s; (c) v = 0.05 m/s; (d) v = 1 m/s.

通過與文獻(xiàn)[22]關(guān)于自由液面的結(jié)果的對(duì)比可以發(fā)現(xiàn), 本研究的混合容器雖然加裝了擋流構(gòu)件, 在不充氣或者低充氣速度時(shí), 可以有效地消除混合過程中可能形成的漩渦的影響, 穩(wěn)定混合空間內(nèi)環(huán)境, 但是當(dāng)充氣速度過高時(shí), 內(nèi)部流場受到強(qiáng)烈擾動(dòng), 自由液面也會(huì)有大幅波動(dòng).

4.2 對(duì)流體速度的影響

為了研究不同的充氣狀態(tài)對(duì)混合空間內(nèi)流體速度的影響, 選取了t= 5 s時(shí), 混合容器內(nèi)徑向位置r= 60 mm, 軸向高度從0到200 mm的速度分布情況, 如圖7所示.從圖7可以看出, 最大速度出現(xiàn)在靠近葉輪的區(qū)域, 底部充氣會(huì)使得出現(xiàn)最大速度的高度上移, 但并不是充氣速度越大最大速度出現(xiàn)的位置越靠上, 四種流體速度分布出現(xiàn)最大速度的高度由低到高依次是充氣速度為v= 0,v= 0.05,v= 0.01和v= 1 m/s時(shí)的混合空間.上述現(xiàn)象主要原因?yàn)楫?dāng)軸向的充氣速度介入流場時(shí), 軸向動(dòng)能和流場的切向動(dòng)能發(fā)生了對(duì)沖能量耗散, 雖然流場的湍流隨機(jī)性增加, 但流動(dòng)的無序性對(duì)整個(gè)流場的總速度有所影響.

圖7 徑向位置r = 60 mm處的軸向高度速度分布 (a) 總速度; (b) 軸向速度; (c) 徑向速度; (d) 切向速度Fig.7.Axial velocity distribution at radial position r = 60 mm: (a) Total velocity; (b) axial velocity; (c) radial velocity; (d) tangential velocity.

從總速度、軸向速度和徑向速度可以看出四種充氣狀態(tài)下的速度分布曲線趨勢走向是相同的, 而對(duì)于切向速度, 在充氣速度為v= 1 m/s的槳葉下方的流體切向速度方向和其他三種情況方向完全相反.顯然, 當(dāng)充氣速度足夠大時(shí), 氣相對(duì)整個(gè)流場的流動(dòng)模式造成較大的擾動(dòng), 剪切流動(dòng)變得復(fù)雜, 增加了流場的湍流混沌特性.此外還可以看出,由于擋流構(gòu)件的存在以及壁面的影響, 可以將切向速度轉(zhuǎn)化為軸向速度和徑向速度, 所以從圖7可以看出切向速度的幅度要小于軸向速度和徑向速度.

4.3 對(duì)顆粒懸浮的影響

圖8、圖9、圖10和圖11是四種充氣狀態(tài)下混合空間內(nèi)三相流的模擬計(jì)算結(jié)果, 截取了部分時(shí)刻的運(yùn)動(dòng)狀態(tài), 這些時(shí)刻基本包含了混合過程中的所有情形, 具有一定的代表性, 其中顆粒的顏色表示顆粒的速度大小.可以很清楚地看到在t= 1 s時(shí), 底部沉積的顆粒在槳葉旋轉(zhuǎn)帶動(dòng)流體的作用下被卷吸起來.初始狀態(tài)時(shí), 攪拌擾動(dòng)流場, 增加了流場的切向流動(dòng)和葉輪底部的軸向上升流運(yùn)動(dòng), 顆粒以較小的速度向上升起.

圖8 v = 0 m/s時(shí)不同時(shí)刻的三相混合系統(tǒng)模擬結(jié)果Fig.8.Simulation results of three-phase stirred system at different time when v = 0 m/s.

圖9 v = 0.01 m/s時(shí)不同時(shí)刻的三相混合系統(tǒng)模擬結(jié)果Fig.9.Simulation results of three-phase stirred system at different time when v = 0.01 m/s.

圖10 v = 0.05 m/s時(shí)不同時(shí)刻的三相混合系統(tǒng)模擬結(jié)果Fig.10.Simulation results of three-phase stirred system at different time when v = 0.05 m/s.

圖11 v = 1 m/s時(shí)不同時(shí)刻的三相混合系統(tǒng)模擬結(jié)果Fig.11.Simulation results of three-phase stirred system at different time when v = 1 m/s.

在t= 1.3 s時(shí), 顆粒群到達(dá)葉輪, 受到高速旋轉(zhuǎn)的作用, 被槳葉打散高速向周圍擴(kuò)散, 同時(shí)可看出, 顆粒群在v= 0.01, 0.05和1 m/s要先于v= 0情況分散開.這是因?yàn)榈撞看禋庠诔跏紶顟B(tài)就增大了流場的軸向剪切流動(dòng)能量, 誘導(dǎo)流體的軸向速度增加,v= 1 m/s時(shí)最為明顯, 如圖7(b)所示.t=1.5 s時(shí), 顆粒受到流場離心力的作用下擴(kuò)散到容器壁面和擋流構(gòu)件, 在兩者的作用下, 顆粒流動(dòng)轉(zhuǎn)變?yōu)樯舷聝蓚€(gè)方向的運(yùn)動(dòng)狀態(tài), 表明當(dāng)旋轉(zhuǎn)流場與擋流構(gòu)件接觸過程時(shí), 流場以切向?yàn)橹鞯牧鲃?dòng)模式遇阻, 切向速度流動(dòng)轉(zhuǎn)變?yōu)樯舷录羟辛鲃?dòng)模式.剩余流場的湍動(dòng)能推動(dòng)部分顆粒作上升流運(yùn)動(dòng), 另一部分顆粒在重力作用向下流動(dòng), 壁面附近的湍動(dòng)能難以對(duì)這部分顆粒提供足夠的驅(qū)動(dòng)動(dòng)能.

隨著時(shí)間的推移, 在t= 1.7 s時(shí), 顆粒到達(dá)自由液面附近, 受到液面振蕩的作用, 分散至容器的各個(gè)位置; 隨后的t= 2和3 s時(shí), 除v= 1 m/s的混合空間外, 其他變化已經(jīng)不明顯, 趨于穩(wěn)態(tài), 可以看出, 速度最大的顆粒分布在葉輪附近, 具有較高速度的顆粒分布在容器的下部, 而越靠近液面處, 顆粒速度越低, 到達(dá)液面時(shí)速度幾乎為0, 此區(qū)域的顆粒群對(duì)氣-液交界面的沖擊非常微小, 顆粒停止上升, 只有擋流構(gòu)件附近的顆粒群受到局部流體漩渦的裹挾作用, 會(huì)繼續(xù)上升, 形成局部的凸液面.此外, 上下兩股流動(dòng)實(shí)現(xiàn)了循環(huán), 對(duì)顆粒分散效果有良好的促進(jìn)作用, 且除v= 1 m/s外, 其他都相對(duì)穩(wěn)定, 顆粒分布都是對(duì)稱的.

從整個(gè)時(shí)間段上三相流的模擬計(jì)算結(jié)果對(duì)比,還可以看出底部充氣使得下部傘狀顆粒群傘柄處要比未充氣狀態(tài)的細(xì), 而v= 1 m/s工況下, 隨著混合的進(jìn)行顆粒受到的作用不均勻, 底部傘狀顆粒群會(huì)被破壞, 運(yùn)動(dòng)狀態(tài)也更加復(fù)雜, 液面波動(dòng)也會(huì)更加劇烈, 難以實(shí)現(xiàn)穩(wěn)態(tài).

為了更準(zhǔn)確、更直觀地描述顆粒分布的均勻性, 引入相對(duì)標(biāo)準(zhǔn)偏差(relative standard deviation, RSD)表征顆粒在液相中的分散效果[36], 在這項(xiàng)工作中, 將液體覆蓋的區(qū)域分為12個(gè)部分(3 ×1 × 4), 即12個(gè)等體積的采樣空間.通過多個(gè)樣本空間內(nèi)顆粒數(shù)目的相對(duì)標(biāo)準(zhǔn)偏差隨時(shí)間的變化作為評(píng)價(jià)指標(biāo).其評(píng)價(jià)公式如下:

其中Xi是采樣空間i中的顆粒數(shù)量,Xavg是Xi的平均值,n為劃分的采樣空間的數(shù)量.從上面的相對(duì)標(biāo)準(zhǔn)偏差RSD評(píng)價(jià)公式可以看出, RSD越小代表顆粒在計(jì)算區(qū)域內(nèi)的顆粒懸浮效果越好.圖12繪制了底部不同充氣工作條件下的RSD隨時(shí)間的變化曲線.

圖12 不同充氣工作條件下的RSD隨時(shí)間的變化Fig.12.RSD changes with time under different aeration conditions.

從圖12可以看出, 在混合容器剛開始工作的一段時(shí)間內(nèi)由于顆粒沉積在底部, 隨時(shí)間移動(dòng)較緩慢, RSD的數(shù)值較大, 顆粒分布極不均勻, 所以RSD曲線有一段上升趨勢, 這段時(shí)間里底部充氣的混合空間由于氣體的作用使得一部分顆粒上升較快, 相對(duì)未充氣的顆粒率先上浮, 所以RSD值要低一些.隨著混合過程的持續(xù)進(jìn)行, 顆粒到達(dá)葉輪高度, 在葉輪的作用下移動(dòng)速度加快, 逐漸分散到混合空間的各個(gè)區(qū)域, RSD曲線呈現(xiàn)快速下降趨勢, 然后趨于小范圍的來回振蕩達(dá)到準(zhǔn)穩(wěn)態(tài)狀態(tài).從圖12也可以發(fā)現(xiàn), 充氣速度為v= 0, 0.01和0.05 m/s的混合空間RSD曲線趨于平緩且處于較低的值, 代表了顆粒分布狀態(tài)趨于穩(wěn)定, 其中v= 0.05 m/s時(shí), 顆粒分散效果最好,v= 0和0.01 m/s次之; 但是充氣速度v= 1 m/s的混合空間的RSD曲線振蕩幅度較大, 因?yàn)槌錃馑俣冗^大導(dǎo)致混合空間內(nèi)局部的不穩(wěn)定性增加, 一定程度上影響了顆粒的懸浮效果.綜上所述, 可以得出選擇適當(dāng)?shù)某錃馑俣瓤梢蕴嵘萜鲀?nèi)的顆粒懸浮效果,提高顆粒分布的均勻性; 而如果充氣速度過大則會(huì)導(dǎo)致混合系統(tǒng)的不穩(wěn)定性, 不利于顆粒懸浮.

5 結(jié) 論

研究具有強(qiáng)剪切過程的氣-液-固三相流混合過程機(jī)理與流場分布, 具有重要的科研價(jià)值和工程應(yīng)用意義.針對(duì)上述問題, 提出了一種研究氣-液-固三相流混合建模與求解方法, 基于VOF-DEM模型建立了流體和顆粒的動(dòng)力學(xué)模型, 進(jìn)行求解, 并將其應(yīng)用在了三相流的混合過程數(shù)值模擬中, 通過結(jié)果分析, 總結(jié)得到了以下結(jié)論:

1) 對(duì)自由液面的變化分析表明, 以較小的速度向帶擋流構(gòu)件的混合空間內(nèi)部充氣可以降低液面的波動(dòng)幅度, 反之, 如果充氣速度過大則會(huì)導(dǎo)致液面波動(dòng)加劇;

2) 對(duì)流體的速度分析表明, 底部充氣會(huì)使得軸向位置上出現(xiàn)最大速度的位置上移, 但并不是充氣速度越大, 該位置就越高, 此外擋流構(gòu)件和壁面的存在可以將流體的切向速度轉(zhuǎn)化為軸向和徑向的速度;

3) 就顆粒懸浮而言, 選擇適當(dāng)?shù)某錃馑俣瓤梢蕴嵘旌峡臻g內(nèi)的顆粒懸浮效果, 有利于相間傳質(zhì), 提高顆粒分布的均勻性, 本次模擬四種狀態(tài)下,v= 0.05 m/s時(shí)顆粒懸浮效果最佳, 而如果充氣速度過大則會(huì)導(dǎo)致混合系統(tǒng)的不穩(wěn)定性, 不利于顆粒懸浮.

基于VOF-DEM模型的三相流混合系統(tǒng)具有高度的復(fù)雜性, 本文在此方面進(jìn)行了有益的嘗試.在理論方面可為VOF-DEM模型的耦合、多相流系統(tǒng)的建模與求解方法等方面研究提供參考; 在技術(shù)方面, 可為工業(yè)混合過程等領(lǐng)域提供更多的技術(shù)解決方案, 具有較好的工程應(yīng)用前景.