付 博，袁希鋼，陳淑勇，劉伯潭，余國琮

氣液界面?zhèn)髻|(zhì)過程Rayleigh對流模擬的格子Boltzmann方法

付 博，袁希鋼，陳淑勇，劉伯潭，余國琮

(天津大學(xué)化工學(xué)院化學(xué)工程聯(lián)合國家重點實驗室，天津 300072)

濃度梯度導(dǎo)致的Rayleigh對流是一種通常在傳質(zhì)過程中觀察到的界面現(xiàn)象．為了模擬乙醇吸收CO2氣液傳質(zhì)過程中界面?zhèn)髻|(zhì)引發(fā)的Rayleigh對流現(xiàn)象，提出了描述界面擾動的隨機(jī)擾動模型，并建立了帶有雙分布模型的二維格子Boltzmann方法(LBM)．模擬二維液相Rayleigh對流過程中，假設(shè)界面上CO2濃度為常數(shù)．通過模擬研究，確定了隨機(jī)擾動模型的2個參數(shù)：局部擾動概率P和濃度擾動大小ρD．模擬結(jié)果表明，當(dāng)0＜ρD≤10-9,kg/m3、10-6≤P≤10-1時，平均液相傳質(zhì)系數(shù)存在一個穩(wěn)定值：(1.09±0.02)×10-5,m/s．通過考察濃度分布結(jié)構(gòu)，分析了Rayleigh對流的時空演化過程．根據(jù)模擬結(jié)果定義并計算了傳質(zhì)增強因子，其證明了Rayleigh對流能夠有效地強化界面?zhèn)髻|(zhì)．

Rayleigh對流；格子Boltzmann方法；界面?zhèn)髻|(zhì)；液相傳質(zhì)系數(shù)；傳質(zhì)增強因子

在相際傳質(zhì)過程中，由物質(zhì)在相界面上的傳遞而導(dǎo)致的濃度梯度會產(chǎn)生密度梯度，若密度梯度與重力作用方向相反，則稱之為界面Rayleigh不穩(wěn)定．此時若遇擾動，則會出現(xiàn)界面失穩(wěn)進(jìn)而引發(fā)的界面流體流動現(xiàn)象，稱之為Rayleigh對流．Rayleigh對流的發(fā)生對界面?zhèn)髻|(zhì)有顯著影響[1]，因而，研究Rayleigh對流對精餾、吸收等化工過程具有重要的意義．以往的實驗研究和理論研究主要集中在采用光學(xué)測試手段如投影法、紋影法、激光全息干涉法等觀察Rayleigh對流結(jié)構(gòu)[2-4]、Rayleigh對流開始的不穩(wěn)定性分析[5-6]、傳質(zhì)速率的測量[7]以及傳質(zhì)系數(shù)關(guān)聯(lián)式的建立[8]．其中實驗研究只能得到傳質(zhì)過程的宏觀信息，而不能得到界面?zhèn)髻|(zhì)的微觀信息；理論研究雖然能夠預(yù)測傳質(zhì)過程中Rayleigh對流發(fā)生時的臨界點，但不能給出Rayleigh對流發(fā)生后的具體傳質(zhì)信息．因此，有必要對Rayleigh對流發(fā)生后詳細(xì)的流動及傳質(zhì)情況進(jìn)行研究．隨著計算機(jī)技術(shù)的發(fā)展，數(shù)值模擬為更好地考察Rayleigh對流及其在強化界面?zhèn)髻|(zhì)中的作用提供了更多可能性，研究結(jié)果將有助于了解Rayleigh對流強化界面?zhèn)髻|(zhì)的機(jī)理．格子Boltzmann方法(lattice Boltzmann method，LBM)作為一種新興的數(shù)值計算方法，由于其在介觀尺度上具有物理條件清晰、邊界處理容易實現(xiàn)、計算效率和數(shù)值精度高等優(yōu)勢，已經(jīng)被成功地應(yīng)用于Rayleigh對流的數(shù)值模擬中．在氣液傳質(zhì)過程中，以往針對液相Rayleigh對流的格子Boltzmann模擬通常假設(shè)界面存在一個或多個始終存在、位置固定并且濃度恒定的高濃度區(qū)域[9-10]．但是，在真實界面?zhèn)髻|(zhì)過程中，界面Rayleigh不穩(wěn)定失穩(wěn)引發(fā)的Rayleigh對流是復(fù)雜的，界面會出現(xiàn)來自各個方面的隨機(jī)擾動，如局部濃度不均勻、物理震動甚至界面的分子熱運動等．因此有必要提出一個求解真實界面?zhèn)髻|(zhì)過程中Rayleigh對流的格子Boltzmann方法．

筆者首先引入一種隨機(jī)擾動模型來描述界面濃度擾動過程，以乙醇吸收CO2界面?zhèn)髻|(zhì)為對象，采用帶有雙分布模型的格子Boltzmann方法模擬液相中二維非穩(wěn)態(tài)Rayleigh對流過程．通過模擬，分析和計算了不同模型參數(shù)下Rayleigh對流的特征：液相傳質(zhì)系數(shù)和濃度分布，通過傳質(zhì)增強因子考察了Rayleigh對流對界面?zhèn)髻|(zhì)的影響．

1 用于Rayleigh對流模擬的格子Boltzmann方法

1.1 雙分布模型

本文采用雙分布模型[11]建立了界面?zhèn)髻|(zhì)引發(fā)Rayleigh對流現(xiàn)象模擬的格子Boltzmann方法．該模型假設(shè)雙組分流體由溶劑A和溶質(zhì)B組成，且B的含量遠(yuǎn)小于A．雙分布模型包含2個粒子分布函數(shù)fAi和gBi，分別用于求解流場和濃度場．無外力項時，基于如圖1所示D2Q9格子模型[12]，則溶劑和溶質(zhì)的LBGK演化方程[13]為

式中：ci為離散速度，m/s；i為基于D2Q9模型的粒子運動方向，i=1，··，9；Δt為單位時間步長，s；τ為無因次松弛時間；fAi(x，t)和gBi(x，t)分別表示t時刻在x位置沿著i方向的溶劑A和溶質(zhì)B的粒子分布函數(shù)．

圖1 二維九速模型(D2Q9模型)Fig.1 Two-dimensional nine-velocity model(D2Q9,model)

溶劑A和溶質(zhì)B的平衡分布函數(shù)[11]為

式中：c為格子速度，c=Δx/Δt，一般取為1；Δx為單位網(wǎng)格步長，m；uA為宏觀速度，m/s；ρ′A為宏觀密度，kg/m3；ρB為溶液中溶質(zhì)的濃度，kg/m3；ωi為權(quán)重系數(shù)，ω1=ω2=ω3=ω4=1/9，ω5=ω6=ω7=ω8= 1/36，ω9=4/9．

通過粒子分布函數(shù)可以計算流場的宏觀變量(密度、速度和濃度)，即

液體的黏度和溶質(zhì)在溶劑中的擴(kuò)散系數(shù)表達(dá)式[11]為

式中：ν為液體的運動黏度，m2/s；D為溶質(zhì)在溶劑中的擴(kuò)散系數(shù)，m2/s．

1.2 外力處理

本文引入Boussinesq近似[14]，即假設(shè)除密度外流場的其他物性(擴(kuò)散系數(shù)、黏度等)為常數(shù)，因而體系所受的外力僅考慮密度變化導(dǎo)致的體積力．本文假設(shè)流體密度是濃度的線性函數(shù)，因此傳質(zhì)導(dǎo)致的濃度的變化是流體所受外部力(即重力)的唯一來源．類似于溫度梯度導(dǎo)致的Rayleigh-Bénard對流中體積力處理方法[14]，則Rayleigh對流中體積力可近似表示為

式中：F為Rayleigh對流中單位體積的體積力，N/m3；g為重力加速度，m/s2；ρ0為溶液的初始濃度，kg/m3；Δρ′為純?nèi)軇┡c飽和溶液之間的密度差，kg/m3；ΔρB為純?nèi)苜|(zhì)與飽和溶液之間的濃度差，kg/m3．

通過修正平衡分布函數(shù)中的變量速度uA來體現(xiàn)外力F對流體的影響．則平衡分布函數(shù)中的速度和流場速度將分別被修正為u*和u[15]，即

1.3 邊界處理

對于流場和濃度場，左右邊界均采用周期邊界；下邊界均采用反彈邊界．對于上邊界，濃度場采用恒濃度邊界[16]；流場采用鏡面對稱邊界[17]．其中鏡面對稱邊界用于處理無變形氣液自由界面時能夠保證水平速度沿豎直方向上的梯度為0，考慮流場被均勻網(wǎng)格覆蓋，其具體做法是假設(shè)上邊界外還有一層虛擬結(jié)點，令虛擬結(jié)點的分布函數(shù)與上邊界結(jié)點的分布函數(shù)沿二者的中心線對稱相等．

1.4 隨機(jī)擾動模型

在界面Rayleigh不穩(wěn)定條件下，一旦有界面擾動，界面會產(chǎn)生Rayleigh對流．界面擾動可能是任何隨機(jī)事件，如機(jī)械振動、分子的熱運動導(dǎo)致的局部濃度擾動等，因而具有隨機(jī)特征．對于在靜止氣液接觸系統(tǒng)的等溫傳質(zhì)過程，可以認(rèn)為擾動僅僅為濃度擾動．因此這里提出了基于濃度擾動的隨機(jī)擾動模型來模型化界面擾動過程，該模型采用2個參數(shù)描述其隨機(jī)特征：局部擾動概率P表示液相界面上每個結(jié)點出現(xiàn)濃度擾動的概率；因擾動導(dǎo)致的濃度凈增量ρD，即濃度擾動大小(設(shè)ρD遠(yuǎn)小于界面濃度ρI)．在LBM模擬中液相界面上每個結(jié)點的濃度0Ιρ表達(dá)式為

式中Βρδ為濃度擾動增量．計算中每個結(jié)點均產(chǎn)生一個隨機(jī)數(shù)，如果該隨機(jī)數(shù)小于P，則Βρδ=±ρD，否則δρΒ=0，且δρΒ的時均值為0．

2 Rayleigh對流的LBM模擬與分析

模擬對象為乙醇溶液吸收CO2過程中，由于CO2通過界面向液體乙醇中擴(kuò)散所導(dǎo)致的Rayleigh對流現(xiàn)象．二維計算區(qū)域如圖2所示，為靜止的液體乙醇，其高為H、寬為L，傳質(zhì)開始前液相為純乙醇．計算區(qū)域的上邊界為氣液自由界面，CO2通過上邊界進(jìn)入液相，下邊界為可視為剛性邊壁．乙醇吸收二氧化碳的氣液傳質(zhì)過程是Rayleigh不穩(wěn)定的和Marangoni穩(wěn)定的[4]．

圖2 Rayleigh對流計算區(qū)域示意Fig.2 Computational domain of Rayleigh convection

模擬基于以下2個假設(shè)：①液膜控制傳質(zhì)過程，界面上氣液兩相達(dá)到平衡狀態(tài)，即界面濃度ρI可視為恒定值[18]，取為4.70,kg/m3；②忽略氣液傳質(zhì)過程中的熱效應(yīng)，且自由界面沒有變形．模擬初始條件：液相主體中各結(jié)點速度均為零，CO2濃度為零．取單位網(wǎng)格步長為5×10-5,m，此時可以獲得可接受的準(zhǔn)確結(jié)果．其中在298.15,K、101.325,kPa條件下，純乙醇及其飽和CO2溶液物性[4,19-20]見表1．表中：0ρ′為純乙醇密度；Bρ?為界面平衡濃度；μ0為純乙醇的動力黏度.

表1 在298.15,K、101.325,kPa條件下純乙醇及其飽和二氧化碳溶液物性Tab.1 Physical properties of ethanol and CO2saturated solution at 298.15,K and 101.325,kPa

2.1 模型參數(shù)的確定

傳質(zhì)系數(shù)是反映傳質(zhì)能力的重要基礎(chǔ)數(shù)據(jù)，可以通過考察模型參數(shù)對傳質(zhì)系數(shù)的影響來確定P和ρD的合適取值區(qū)間．取計算區(qū)域高H=5,mm、寬L= 5,mm，網(wǎng)格劃分為100×100．

t時間內(nèi)平均液相傳質(zhì)系數(shù)kavg,t(m/s)可表示為

式中：V為溶液體積，m3；AI為界面面積，m2；ρt為t時刻溶液中溶質(zhì)的平均濃度，kg/m3．

圖3給出了不同局部擾動概率P時t=300,s內(nèi)平均液相傳質(zhì)系數(shù)kavg隨濃度擾動大小Dρ的變化關(guān)系．

圖3 不同P時kavg隨ρD的變化關(guān)系Fig.3 Variations of kavgwith ρDfor different P

從圖3中可以看出，給定P，kavg隨著ρD的減小呈現(xiàn)出先減小后趨于穩(wěn)定的變化趨勢：當(dāng)ρD＞10-9kg/m3時，kavg隨著ρD的減小而減小；當(dāng)0＜ρD≤10-9kg/m3時，kavg與ρD的取值基本無關(guān)．在0＜ρD≤10-9kg/m3的條件下，當(dāng)10-6≤P≤10-1時，kavg存在一個穩(wěn)定值，為(1.09±0.02)×10-5m/s，則相應(yīng)的平均液相傳質(zhì)通量為(5.12±0.094)×10-5kg/(m2·s)，與文獻(xiàn)[21]實驗測量結(jié)果(300,s內(nèi)平均液相傳質(zhì)通量約為6×10-5kg/(m2·s))吻合較好，從而說明了模型參數(shù)P和ρD取值范圍的合理性．

2.2 Rayleigh對流結(jié)構(gòu)

取計算區(qū)域高H=10,mm、寬L=20,mm，網(wǎng)格劃分為200×400．圖4給出了50,s、55,s、60,s和65,s時不同模型參數(shù)下溶質(zhì)的濃度分布．當(dāng)P=0.1、ρD= 10-12,kg/m3時，從圖4(a)～(d)可以觀察到Rayleigh對流的時空演化過程：在有界面擾動條件下，溶質(zhì)首先在近界面發(fā)生擴(kuò)散，約50,s左右，界面由于Rayleigh不穩(wěn)定開始發(fā)生變形(見圖4(a))；55,s時界面對流已經(jīng)在多個Rayleigh不穩(wěn)定區(qū)域開始發(fā)生(見圖4(b))；隨后對流發(fā)展成羽狀對流結(jié)構(gòu)，界面處高濃度液體以羽狀對流結(jié)構(gòu)逐步向下發(fā)展至液相主體(見圖4(c)和(d))．這種Rayleigh對流的時空演化特征同實驗觀察結(jié)果基本一致[2-4]．

圖4 50,s、55,s、60,s和65,s時不同模型參數(shù)下溶質(zhì)的濃度分布Fig.4 Transient contours of solute concentration for different model parameters after 50,s，55,s，60,s，65,s

對比圖4(d)和圖4(e)可知，當(dāng)局部擾動概率P一定時，羽狀對流結(jié)構(gòu)的數(shù)量、大小及產(chǎn)生位置與濃度擾動大小ρD基本無關(guān)．對比圖4(e)和圖4(f)可知，當(dāng)ρD一定時，P對羽狀對流結(jié)構(gòu)的產(chǎn)生位置有一定的影響，但對其數(shù)量及大小基本無影響．

2.3 傳質(zhì)增強因子

發(fā)生Rayleigh對流時，液相傳質(zhì)系數(shù)比只依靠擴(kuò)散傳質(zhì)時大，因此引出傳質(zhì)增強因子這一概念，將其定義為存在Rayleigh對流時的液相傳質(zhì)系數(shù)與只有擴(kuò)散傳質(zhì)存在時(由滲透理論計算)的液相傳質(zhì)系數(shù)之比[8]．

根據(jù)模擬獲得的濃度場，t時刻瞬時液相傳質(zhì)系數(shù)kins,t(m/s)計算式為

式中ΔT為時間間隔，取值為0.1,s．

由滲透理論預(yù)測的瞬時液相傳質(zhì)系數(shù)k0(m/s)ins,t可表示為

則瞬時傳質(zhì)增強因子可表示為

取計算區(qū)域高H=10,mm、寬L=20,mm，網(wǎng)格劃分為200×400．圖5給出了P=0.1、ρD=10-12,kg/m3條件下瞬時傳質(zhì)增強因子Fins隨時間的變化關(guān)系．由圖5可知，瞬時傳質(zhì)增強因子介于1～7之間，因此Rayleigh對流可以大大提高氣液傳質(zhì)速率．

圖5 P=0.1、ρD=10-12,kg/m3條件下Fins隨時間t變化關(guān)系Fig.5 Variation of Finswith time t at P=0.1,and ρD=10-12,kg/m3

3 結(jié) 語

通過引入隨機(jī)擾動模型，建立了模擬Rayleigh對流的帶有雙分布模型的二維格子Boltzmann方法．運用該方法對靜止乙醇吸收CO2傳質(zhì)過程中液相Rayleigh對流現(xiàn)象進(jìn)行模擬，結(jié)果發(fā)現(xiàn)，300,s內(nèi)平均液相傳質(zhì)系數(shù)kavg隨著ρD的減小呈現(xiàn)出先減小后趨于穩(wěn)定的變化趨勢；當(dāng)0＜ρD≤10-9,kg/m3、10-6≤P≤10-1時，kavg存在一個穩(wěn)定值(1.09± 0.02)×10-5,m/s．通過考察濃度分布結(jié)構(gòu)，分析了Rayleigh對流的時空演化過程．另外，傳質(zhì)增強因子的模擬結(jié)果證明了Rayleigh對流能夠有效地提高傳質(zhì)速率．因此，建立的LBM方法可有效模擬界面?zhèn)髻|(zhì)引發(fā)的Rayleigh對流現(xiàn)象，有助于加強對真實氣液傳質(zhì)過程中Rayleigh對流的理解，為進(jìn)一步探討強化界面?zhèn)髻|(zhì)的機(jī)理和途徑提供幫助．

［1］ Grahn A. Two-dimensional numerical simulations of Marangoni-Bénard instabilities during liquid-liquid mass transfer in a vertical gap[J]. Chemical Engineering Science，2006，61(11)：3586-3592.

［2］ 沙 勇，李樟云，林芬芬，等. 氣液傳質(zhì)界面湍動現(xiàn)象投影觀察[J]. 化工學(xué)報，2010，61(4)：844-847. Sha Yong，Li Zhangyun，Lin Fenfen，et al. Shadowgraph observation on interfacial turbulence phenomena in gas-liquid mass transfer[J]. CIESC Journal，2010，61(4)：844-847(in Chinese).

［3］ Barbosa J R，Ortolan M A. Experimental and theoretical analysis of refrigerant absorption in lubricant oil[J]. HVAC and R Research，2008，14(1)：141-158.

［4］ Okhotsimskiis A，Hozawa M. Schlieren visualization of natural convection in binary gas-liquid systems[J]. Chemical Engineering Science，1998，53(14)：2547-2573.

［5］ Sun Z F. Onset of Rayleigh-Bénard-Marangoni convection with time-dependent nonlinear concentration profiles[J]. Chemical Engineering Science，2012，68 (1)：579-594.

［6］ Tan K K，Tey B T，Tan Y W. Onset of natural convection in gas-gas system induced by bottom-up transient mass diffusion[J]. Engineering Applications of Computational Fluid Mechanics，2010，4(4)：475-482.

［7］ 陳 煒. 氣液界面Rayleigh-Bénard-Marangoni對流現(xiàn)象實驗測量及傳質(zhì)研究[D]. 天津：天津大學(xué)化工學(xué)院，2010. Chen Wei. Experimental Measurement of Gas-Liquid Interfacial Rayleigh-Bénard-Marangoni Convection and Mass Transfer[D]. Tianjin：School of Chemical Engineering，Tianjin University，2010 (in Chinese).

［8］ Sun Z F，Yu K T，Wang S Y，et al. Absorption anddesorption of carbon dioxide into and from organic solvents：Effect of Rayleigh and Marangoni instability[J]. Ind Eng Chem Res，2002，41(7)：1905-1913.

［9］ 付 博，袁希鋼，陳淑勇，等. Rayleigh對流及其對界面?zhèn)髻|(zhì)影響模擬的格子Boltzmann方法[J]. 化工學(xué)報，2011，62(11)：2993-3000. Fu Bo，Yuan Xigang，Chen Shuyong，et al. Rayleigh convection and its effect on interfacial mass transfer by lattice Boltzmann simulation[J]. CIESC Journal，2011，62(11)：2993-3000(in Chinese).

［10］ Fu Bo，Yuan Xigang，Liu Botan，et al. Characterization of Rayleigh convection in interfacial mass transfer by lattice Boltzmann simulation and experimental verification[J]. Chinese Journal of Chemical Engineering，2011，19(5)：845-854.

［11］ Inamuro T，Yoshino M，Inoue H，et al. A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat-transfer problem[J]. Journal of Computational Physics，2002，179(1)：201-215.

［12］ Qian Y H，D'Humieres D，Lallemand P. Lattice BGK models for Navier-Stokes equation[J]. Europhysics Letters，1992，17(6)：479-484.

［13］ Bhatnagar P L，Gross E P，Krook M. A model for collision processes in gases. Ⅰ. small amplitude processes in charged and neutral one-component systems[J]. Physical Review，1954，94(3)：511-525.

［14］ Shan X W. Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method[J]. Physical Review E，1997，55(3)：2780-2788.

［15］ Buick J M，Greated C A. Gravity in a lattice Boltzmann model[J]. Physical Review E，2000，61(5)：5307-5320.

［16］ Sukop M C，Thorne D T Jr. Lattice Boltzmann Modeling：An Introduction for Geoscientists and Engineers[M]. Netherlands：Springer，2006.

［17］ Shi Y，Zhao T S，Guo Z L. Finite difference-based lattice Boltzmann simulation of natural convection heat transfer in a horizontal concentric annulus[J]. Computers and Fluids，2006，35：1-15.

［18］ Tan K K，Thorpe R B. The onset of convection induced by buoyancy during gas diffusion in deep fluids[J]. Chemical Engieering Science，1999，54：4179-4187.

［19］ Arce A，Arce A Jr，Rodil E，et al. Density，refractive index，and speed of sound for 2-ethoxy-2-methylbutane+ ethanol+water at 298. 15,K[J]. Journal of Chemical and Engineering Data，2000，45(4)：536-539.

［20］ Takahashi M，Kobayashi Y. Diffusion coefficients and solubilities of carbon dioxide in binary mixed solvents[J]. Journal of Chemical and Engineering Data，1982，27(3)：328-331.

［21］ 林芬芬. 界面湍動現(xiàn)象的實驗研究[D]. 廈門：廈門大學(xué)化學(xué)化工學(xué)院，2009. Lin Fenfen. Experimental Study on the Interfacial Turbulence Phenomena[D]. Xiamen：College of Chemistry and Chemical Engineering，Xiamen University，2009(in Chinese).

Simulation of Rayleigh Convection in Gas-Liquid Interfacial Mass Transfer Process by Lattice Boltzmann Method

FU Bo，YUAN Xi-gang，CHEN Shu-yong，LIU Bo-tan，YU Kuo-tsung
(State Key Laboratory of Chemical Engineering，School of Chemical Engineering and Technology，Tianjin University，Tianjin 300072，China)

Concentration gradient driven Rayleigh convection is an interfacial phenomenon usually observed in mass transfer processes. A random perturbation model characterizing interfacial disturbance was proposed and a twodimensional lattice Boltzmann method(LBM)with a double distribution model was implemented for simulating the interfacial mass transfer induced Rayleigh convection in gas-liquid mass transfer process of CO2absorption into liquid ethanol. In the simulation of the two-dimensional liquid-phase Rayleigh convection，CO2concentration at the interface was assumed to be constant. Through simulation，two parameters for the random perturbation model—the local perturbation possibility P and the perturbation magnitude ρD—were identified. The simulated results show that an average liquid-phase mass transfer coefficient has a stable value of(1.09±0.02)×10-5,m/s with the parameter values of 0＜ρD≤10-9kg/m3and 10-6≤P≤10-1. Through investigating concentration contour patterns，the temporal-spatial evolution of the convection was analyzed. A mass transfer enhancement factor was defined and evaluated based on the simulated results，and the interfacial mass transfer was shown to be effectively enhanced by the Rayleigh convection.

Rayleigh convection；lattice Boltzmann method(LBM)；interfacial mass transfer；liquid-phase mass transfer coefficient；mass transfer enhancement factor

TQ021.4

A

0493-2137(2012)07-0585-06

2012-01-14；

2012-04-19.

國家自然科學(xué)基金資助項目(20736005).

付 博（1984— ），男，博士研究生，fubo@tju.edu.cn.

袁希鋼，yuanxg@tju.edu.cn.