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

(浙江工業(yè)大學機械工程學院, 特種裝備制造與先進加工技術教育部/浙江省重點實驗室, 杭州 310014)

1 引 言

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

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

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

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

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

2 三相流動力學模型

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

2.1 VOF模型

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

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

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

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

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

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

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

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

2.2 DEM模型

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

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

法向接觸力:

切向接觸力:

切向摩擦力矩:

滾動摩擦力矩:

顆粒-流體作用力:

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

2.3 耦合模型

為了獲得流體與顆粒間準確的相互作用力, 提出了一種相間耦合解法—多孔模型來描述精確的顆粒運動軌跡, 其計算表達式如下:

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

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

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

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

3 數(shù)值計算

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

研究所選取的物理混合空間為帶擋流板的半圓底容器, 混合執(zhí)行構件為葉輪, 結構如圖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, 底部進氣口直徑d2= 14 mm, 高度ha= 4 mm, 擋流板高度hb= 11T/10, 寬度為T/10.

圖2 混合空間結構示意圖Fig.2.Diagram of mixed space structure.

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

圖3 網(wǎng)格劃分 (a) 靜子區(qū)域網(wǎng)格; (b) 轉子區(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模型計算,選擇顯式的時間離散格式, 湍流模型使用標準的湍動能-耗散(k-ε)模型, 該模型具有較高的物理可靠性, 可為混合湍流過程提供精確解, 近壁區(qū)域選擇標準壁面函數(shù).混合容器頂部設置為壓力出口邊界條件, 容器壁面為無滑移壁面邊界條件, 底部為充氣管道, 并采用速度入口邊界條件.數(shù)值求解方法使用coupled方案, 該解法耦合了動量、能量及組分方程, 能比較快地得到收斂解, 動量離散格式、湍動能和湍動能耗散率離散格式均采用二階迎風以獲得精確的解[29-31], 體積分數(shù)離散格式使用分段線性界面重構(piecewise linear interface construction, PLIC)算法, 這種方法是精度最高的一種[24,29], 監(jiān)視器收斂殘差均為10—6.

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

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

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

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

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

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

4 結果與討論

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

4.1 對自由液面的影響

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

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

圖5 t = 5 s時四種充氣條件下混合空間內的自由液面 (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.

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

圖6 t = 5 s時, 四種充氣條件下混合空間內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.

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

4.2 對流體速度的影響

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

圖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)下的速度分布曲線趨勢走向是相同的, 而對于切向速度, 在充氣速度為v= 1 m/s的槳葉下方的流體切向速度方向和其他三種情況方向完全相反.顯然, 當充氣速度足夠大時, 氣相對整個流場的流動模式造成較大的擾動, 剪切流動變得復雜, 增加了流場的湍流混沌特性.此外還可以看出,由于擋流構件的存在以及壁面的影響, 可以將切向速度轉化為軸向速度和徑向速度, 所以從圖7可以看出切向速度的幅度要小于軸向速度和徑向速度.

4.3 對顆粒懸浮的影響

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

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

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

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

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

在t= 1.3 s時, 顆粒群到達葉輪, 受到高速旋轉的作用, 被槳葉打散高速向周圍擴散, 同時可看出, 顆粒群在v= 0.01, 0.05和1 m/s要先于v= 0情況分散開.這是因為底部吹氣在初始狀態(tài)就增大了流場的軸向剪切流動能量, 誘導流體的軸向速度增加,v= 1 m/s時最為明顯, 如圖7(b)所示.t=1.5 s時, 顆粒受到流場離心力的作用下擴散到容器壁面和擋流構件, 在兩者的作用下, 顆粒流動轉變?yōu)樯舷聝蓚€方向的運動狀態(tài), 表明當旋轉流場與擋流構件接觸過程時, 流場以切向為主的流動模式遇阻, 切向速度流動轉變?yōu)樯舷录羟辛鲃幽Ｊ?剩余流場的湍動能推動部分顆粒作上升流運動, 另一部分顆粒在重力作用向下流動, 壁面附近的湍動能難以對這部分顆粒提供足夠的驅動動能.

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

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

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

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

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

從圖12可以看出, 在混合容器剛開始工作的一段時間內由于顆粒沉積在底部, 隨時間移動較緩慢, RSD的數(shù)值較大, 顆粒分布極不均勻, 所以RSD曲線有一段上升趨勢, 這段時間里底部充氣的混合空間由于氣體的作用使得一部分顆粒上升較快, 相對未充氣的顆粒率先上浮, 所以RSD值要低一些.隨著混合過程的持續(xù)進行, 顆粒到達葉輪高度, 在葉輪的作用下移動速度加快, 逐漸分散到混合空間的各個區(qū)域, RSD曲線呈現(xià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時, 顆粒分散效果最好,v= 0和0.01 m/s次之; 但是充氣速度v= 1 m/s的混合空間的RSD曲線振蕩幅度較大, 因為充氣速度過大導致混合空間內局部的不穩(wěn)定性增加, 一定程度上影響了顆粒的懸浮效果.綜上所述, 可以得出選擇適當?shù)某錃馑俣瓤梢蕴嵘萜鲀鹊念w粒懸浮效果,提高顆粒分布的均勻性; 而如果充氣速度過大則會導致混合系統(tǒng)的不穩(wěn)定性, 不利于顆粒懸浮.

5 結 論

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

1) 對自由液面的變化分析表明, 以較小的速度向帶擋流構件的混合空間內部充氣可以降低液面的波動幅度, 反之, 如果充氣速度過大則會導致液面波動加劇;

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

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

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