格子Boltzmann方法模擬咸水吸收CO2的Rayleigh對流過程

付博，田建昌，劉菊，張潤葉，陳慕華，朱新寶

（南京林業(yè)大學(xué)化學(xué)工程學(xué)院，江蘇南京210037）

化石燃料燃燒導(dǎo)致CO2大量排放引起的溫室效應(yīng)嚴(yán)重威脅著人類賴以生存的地球環(huán)境，遏制全球氣候變暖、減少大氣中CO2濃度成為全球關(guān)注的熱點。二氧化碳捕集與封存（CO2capture and storage，CCS）是CO2大規(guī)模減排應(yīng)對氣候變化的重要手段，被認(rèn)為是最有潛力和最可行的CO2減排方式[1]。而咸水層在全球分布廣泛、碳封存能力強(qiáng)，受到了廣泛關(guān)注[2]。研究CO2在咸水層的遷移過程，對于設(shè)計咸水層封存項目至關(guān)重要。在封存階段，為CO2溶解和遷移提供動力的是分子擴(kuò)散和界面對流[3-4]。即在氣液相界面處CO2通過分子擴(kuò)散的方式進(jìn)入咸水，同時界面處已溶解CO2的遷移速率制約著CO2進(jìn)一步的溶解。與分子擴(kuò)散相比，界面處的對流傳質(zhì)可明顯加快CO2的溶解速率[5]。在CO2溶解過程中，液體密度明顯增大，而表面張力減小，因此本文中涉及的對流現(xiàn)象可視為只由密度梯度引發(fā)的Rayleigh對流[6]。CO2溶解處密度差持續(xù)積累導(dǎo)致重力不穩(wěn)定，突破黏性阻力發(fā)生沉降，而不含CO2的咸水則會向上流動，由此產(chǎn)生對流，加快了CO2進(jìn)入溶液的速率[7]。在不同地質(zhì)條件下，除了溫度和壓力外，咸水層的鹽濃度差異也會影響CO2溶解和遷移速率，因此需要了解CO2在溶于咸水過程中對流傳質(zhì)的規(guī)律以及鹽濃度對傳質(zhì)速率的影響。

目前對于CO2在咸水中的溶解和對流傳質(zhì)過程的研究，主要有實驗法和數(shù)值模擬法。Khosrokhavar等[8]采用PVT筒檢測了CO2溶于鹽水過程中密度梯度的變化，整體了解CO2在咸水中的溶解特性。為了更直接地觀察和研究自然對流的整個過程，Kneafsey等[9]利用垂直的Hele-Shaw盒子觀測了常壓下CO2指狀Rayleigh對流結(jié)構(gòu)的形成和發(fā)展過程，同時指出，外界微小的擾動都會造成不同Rayleigh對流結(jié)構(gòu)的產(chǎn)生。Soroush等[10]也通過Hele-Shaw 盒子實驗，驗證了Rayleigh對流結(jié)構(gòu)對傳質(zhì)通量的提高具有很強(qiáng)的促進(jìn)作用。Abad等[11]和Khosrokhavar等[12]的實驗結(jié)果表明，與在純水中CO2的溶解過程相比，CO2在咸水中的對流傳質(zhì)過程明顯減弱。隨著計算機(jī)模擬技術(shù)的發(fā)展，通過建立適當(dāng)?shù)臄?shù)學(xué)模型來模擬Rayleigh對流的產(chǎn)生和結(jié)構(gòu)變化，可以更清晰直觀地了解該過程的機(jī)理。MacMinn等[13]和Xu等[14]基于穩(wěn)定性分析的模擬，研究了CO2注入咸水層后從密度梯度增加到垂直對流發(fā)生的過程，并確定了Rayleigh對流開始的臨界時間，但未對之后的混合過程進(jìn)行模擬。Pruess等[15]基于高網(wǎng)格精度數(shù)值模擬，研究了CO2的封存過程中溶解-擴(kuò)散-對流的過程，較好地模擬了Rayleigh對流結(jié)構(gòu)從產(chǎn)生到發(fā)展至下邊界之前的過程，但該模擬方法只考慮了密度差引起的體積力的影響，對于咸水中鹽離子對溶解速率的影響沒有進(jìn)一步探索。

本文在前人研究的基礎(chǔ)上，采用格子Boltzmann方法（LBM）模擬研究CO2在溶解于咸水過程中形成的對流結(jié)構(gòu)特點和傳質(zhì)規(guī)律以及鹽濃度對該過程的影響。與其他數(shù)值模擬方法相比，LBM具有數(shù)值精度和計算效率高、物理條件清晰、邊界條件易處理等優(yōu)點，可以在介觀尺度上有效模擬流體流動和傳質(zhì)過程[16]。本文根據(jù)實際的傳質(zhì)過程特點，建立了包含密度分布函數(shù)和濃度分布函數(shù)的雙分布模型，并引入靜電力和體積力等外力項。模擬了在氣液界面處存在多個飽和CO2擴(kuò)散源時，Rayleigh 對流從產(chǎn)生到發(fā)展再到穩(wěn)定的全過程，并對CO2在純水和不同鹽濃度溶液中的溶解過程進(jìn)行了模擬，驗證了前人理論分析和實驗中所述的鹽濃度對Rayleigh對流過程的抑制作用，模擬結(jié)果可為CO2捕集與封存提供理論指導(dǎo)。

1 建立Rayleigh對流模擬的LBM模型

1.1 速度場

本文采用Inamuro 等[17]提出的雙分布模型，假設(shè)雙組分流體由溶劑A和溶質(zhì)B組成，且B的含量遠(yuǎn)小于A，因此B 對流場變化的貢獻(xiàn)可忽略不計。計算過程中采用D2Q9格子模型[18]，見圖1。

溶劑A的無外力LBGK演化方程[19]見式(1)。

式中，fAi(x,t)為t時刻在x位置上沿i方向溶劑A的粒子分布函數(shù)；τA為A 粒子兩次碰撞的時間間隔，即A 粒子分布函數(shù)趨于平衡分布函數(shù)的松弛時間。

圖1 二維九速模型（D2Q9模型）

溶劑A的平衡分布函數(shù)見式(2)。

式中，c=Δx/Δt 為格子速度，即單位網(wǎng)格步長和單位時間步長的比值；權(quán)重系數(shù)ωi和離散速度ci分別見式(3)和式(4)。

可由粒子分布函數(shù)計算流場的宏觀密度和宏觀速度，分別見式(5)和式(6)。

液體黏度見式(7)。

1.2 濃度場

溶質(zhì)B的粒子分布情況即為該擴(kuò)散過程的濃度場的變化，無外力LBGK演化方程見式(8)。

式中，gBi(x,t)表示t 時刻在x 位置上沿i 方向溶質(zhì)B的粒子分布函數(shù)，即濃度分布函數(shù)；τB是溶質(zhì)B的松弛時間。

溶質(zhì)B的平衡分布函數(shù)見式(9)。

宏觀濃度見式(10)。

溶質(zhì)擴(kuò)散系數(shù)[20]見式(11)。

1.3 雙分布模型

本文采用Buick 等[21]提出的混合模型，通過修改碰撞項和平衡分布函數(shù)，將速度場和濃度場耦合，達(dá)到處理外力的目的。

在該體系中，鹽、水和二氧化碳之間存在3種相互作用力：①離子和水分子間的靜電力；②二氧化碳和水分子之間的范德華力；③離子和二氧化碳分子間的靜電力。Chalbaud 等[22]和Aggelopoulos等[23-24]研究表明，在該體系中，離子與二氧化碳之間的作用力占主導(dǎo)，其他兩種力遠(yuǎn)小于該力。流場中鈉離子、氯離子與CO2分子間的靜電力作用產(chǎn)生的勢能[25]分別見式(12)和式(13)。

式中，z1和z2分別為鈉離子和氯離子的帶電荷數(shù)；e為單位電荷的帶電量，取1.602177×10-19C；α為二氧化碳的極化率；ε0為真空介電常數(shù)；r 為離子與分子中心之間的距離，為計算方便，距離統(tǒng)一取臨界間距0.5nm[26-27]。

除了鹽的靜電力外，咸水吸收CO2會引起體系密度的變化，由此產(chǎn)生的體積力也是Rayleigh對流傳質(zhì)過程中的重要驅(qū)動力之一[28]。本文采用Boussinesq 近似[29]來處理體積力，即假設(shè)流體密度是濃度的線性函數(shù)，由此產(chǎn)生的體積力見式(14)。

因此，單位體積的溶質(zhì)粒子在格子結(jié)點上受到的總作用力勢能見式(15)。

修改后的碰撞方程見式(16)。

式中，Bi為速度系數(shù)，其目的是加入受力影響時，保證質(zhì)量和動量守恒。其中B1=B2=B3=B4=B9=6，B5=B6=B7=B8=24。

修正后平衡狀態(tài)函數(shù)中的速度見式(17)。

外力修正后的流場速度[30]見式(18)。

修正后的平衡狀態(tài)方程分別見式(19)和式(20)。

1.4 邊界處理

計算區(qū)域如圖2 所示，區(qū)域的高為H，寬為L，溫度和壓力設(shè)置為293.15K、101kPa，并忽略傳質(zhì)過程中溫度的變化。初始時刻，液相為靜止的、均勻的NaCl 水溶液，各結(jié)點速度和CO2濃度均為0。區(qū)域下邊界為靜止無滑移壁面，采用標(biāo)準(zhǔn)反彈邊界[31]；為防止左右兩側(cè)對Rayleigh 對流結(jié)構(gòu)的發(fā)展產(chǎn)生影響，左右采用周期邊界[19]；上邊界是無變形氣液自由界面，采用鏡面對稱邊界[32]；擴(kuò)散源處CO2的濃度為飽和濃度，采用恒濃度邊界處理[33]。此時，CO2及其飽和溶液的物性見表1。

圖2 Rayleigh對流計算區(qū)域

表1 不同NaCl濃度咸水飽和二氧化碳溶液物性[34-37]

2 Rayleigh對流的LBM模擬和分析

本文采用Matlab 軟件編程，計算過程主要包括粒子碰撞和遷移[38]：①初始化流場每個結(jié)點的變量及分布函數(shù)；②計算每個結(jié)點的宏觀物理量，代入式(19)和式(20)，由此得到粒子的平衡分布函數(shù)；③在每個結(jié)點上進(jìn)行碰撞計算；④在流場內(nèi)部結(jié)點中，將粒子遷移到相鄰結(jié)點；⑤在邊界結(jié)點上進(jìn)行邊界條件處理；⑥得到經(jīng)過一次碰撞和遷移后的宏觀量，并作為下一次迭代的初始值，經(jīng)過多次迭代計算，最終得到傳質(zhì)過程中的速度場和濃度場。

2.1 單個及多個溶質(zhì)擴(kuò)散源的模擬

2.1.1 單個溶質(zhì)擴(kuò)散源的模擬

對于單個CO2擴(kuò)散源的傳質(zhì)過程，設(shè)定的計算區(qū)域高H=5mm、寬L=5mm，網(wǎng)格劃分為100×100，則單位網(wǎng)格長度為Δx=0.05mm。圖3 給出了單個CO2擴(kuò)散源溶解于1.0mol/L 的NaCl 溶液過程中，溶質(zhì)濃度分布情況的模擬結(jié)果及相關(guān)實驗圖像。

由圖3 的濃度變化可知，在只有單個CO2擴(kuò)散源時，指狀對流結(jié)構(gòu)不受墻壁和其他對流結(jié)構(gòu)的影響，界面處CO2濃度較高的液體沿重力方向流動，逐漸發(fā)展為對稱的Rayleigh 對流結(jié)構(gòu)，頭部圓潤，尾部瘦小，具有明顯的流線型輪廓，這與圖3(c)中Thomas等[39]實驗觀察到的對流結(jié)構(gòu)十分相近。

2.1.2 多個溶質(zhì)擴(kuò)散源的模擬

對于多個擴(kuò)散源的傳質(zhì)過程，設(shè)定的計算區(qū)域高H=10mm、寬L=20mm，網(wǎng)格劃分為100×200，則單位網(wǎng)格長度Δx=0.1mm。n個CO2擴(kuò)散源均勻分布在氣液界面，圖4 給出了n=5 時CO2溶解于1.0mol/L的NaCl溶液過程中不同時刻溶質(zhì)的濃度分布情況。

由圖4的濃度分布可知，指狀對流結(jié)構(gòu)經(jīng)歷了從生長到發(fā)展再到觸底反彈的過程。t=50s左右時，在氣液界面處形成了5個指狀對流結(jié)構(gòu)，并且開始相互誘導(dǎo)。t=100s時，兩側(cè)的指狀對流結(jié)構(gòu)（擴(kuò)散源1和2，4和5）從根部合并，完全融合并形成較大對流結(jié)構(gòu)，繼續(xù)向下和兩側(cè)發(fā)展。t=150s時，原來的5個指狀對流結(jié)構(gòu)在根部合并，對流結(jié)構(gòu)末端向上折返，為大對流結(jié)構(gòu)向下發(fā)展提供條件，這種分支和逆向流動可促進(jìn)傳質(zhì)，與沙勇等[40]通過雙光路紋影儀實驗觀察到的半圓形小渦流結(jié)構(gòu)非常相似。下半部分的兩個分支到達(dá)下邊界后，在壁面的影響下向上返回，并繼續(xù)在橫向發(fā)展。t=200s 時，原先兩個分支向下發(fā)展融合成一體，Rayleigh 對流結(jié)構(gòu)發(fā)展充分，達(dá)到穩(wěn)定狀態(tài)。

2.1.3 壁面阻力的影響

本文模擬計算的對象是類似于Hele-Shaw盒子實驗的設(shè)定區(qū)域，兩塊垂直的薄板間隙非常小，因此可以將實際三維的薄板間的Rayleigh對流過程近似為二維層流流動，服從平面泊肅葉流規(guī)律，受到的壁面阻力可表示為式(21)[41]。

式中，δ 為壁面位置；μ 為溶液的動力黏度。將式(21)代入式(17)和式(18)，分別修正平衡速度和流場速度，由此實現(xiàn)壁面阻力對該過程的影響。模擬結(jié)果如圖5 所示，擴(kuò)散源數(shù)n=5，比較在t=150s時不同壁面位置處x-y截面的濃度分布圖可知，隨著與壁面的靠近，壁面阻力帶來的影響越來越大，Rayleigh對流過程受到抑制。

2.2 鹽濃度對Rayleigh對流結(jié)構(gòu)的影響

在293.15K、101kPa 下，分別在純水和摩爾濃度 為0.5mol/L、 1.0mol/L、 1.5mol/L、 2.0mol/L、2.5mol/L 和3.0mol/L 的溶液中進(jìn)行了二氧化碳在溶液中的自然對流傳質(zhì)模擬。設(shè)定的計算區(qū)域高H=10mm、寬L=20mm，網(wǎng)格劃分為100×200，單位網(wǎng)格長度Δx=0.1mm。圖6 給出了n=7 時不同鹽濃度和不同時刻溶質(zhì)的濃度分布情況。

圖5 n=5、t=150s時不同壁面位置處x-y截面的濃度分布

觀察圖6 中的濃度變化可知，在t=50s 時，純水以及鹽濃度c0≤1.0mol/L 的溶液中，氣液界面附近出現(xiàn)了明顯的指狀對流結(jié)構(gòu)，并且結(jié)構(gòu)之間相互誘導(dǎo)，開始發(fā)生形變；而在鹽濃度c0＞1.0mol/L的溶液中，指狀對流結(jié)構(gòu)發(fā)展緩慢，結(jié)構(gòu)之間幾乎沒有影響。在t=100s 時，純水以及鹽濃度c0=0.5mol/L 的溶液中，小指狀對流結(jié)構(gòu)在根部合并，合并后大指狀對流結(jié)構(gòu)迅速向下發(fā)展并往橫向延伸，未參與合并的指狀對流結(jié)構(gòu)發(fā)展緩慢；在鹽濃度c0＞0.5mol/L 的溶液中，小指狀對流結(jié)構(gòu)還未發(fā)生融合。在t=200s 時，所有溶液中的對流結(jié)構(gòu)在近氣液界面處完全融合，向下發(fā)展的大指狀對流結(jié)構(gòu)觸底反彈。其中，純水以及鹽濃度c0≤1.0mol/L 的溶液中，大指狀對流結(jié)構(gòu)已發(fā)展至底部并趨于融合，致使底部平均濃度較高，以向上發(fā)展為主；而在鹽濃度c0＞1.0mol/L 的溶液中，部分指狀對流結(jié)構(gòu)還未發(fā)展至溶液底部。

圖6 的模擬結(jié)果表明，隨著溶液鹽濃度的增大，指狀對流結(jié)構(gòu)出現(xiàn)的時間會推遲，并使后續(xù)的對流傳質(zhì)過程減弱。這與楊盼瑞等[42]通過實驗觀測到的二氧化碳在咸水中鹽濃度對擴(kuò)散系數(shù)的影響情況一致。由此可知，對于不同鹽濃度的地質(zhì)條件，CO2封存過程必須充分考慮鹽濃度對CO2擴(kuò)散和對流傳質(zhì)的影響。

2.3 模擬結(jié)果分析

2.3.1 對流發(fā)生時間及流場速度分析

在1.0mol/L 的NaCl 溶液中，不同擴(kuò)散源個數(shù)時對流的開始時間模擬結(jié)果如圖7 所示。由圖可知，隨著擴(kuò)散源個數(shù)的增加，Rayleigh 對流發(fā)生的時間越來越長。至n=40 時，對流發(fā)生時間tonset=52.4s，與Loodts等[43]對連續(xù)擴(kuò)散源模型的對流發(fā)生時間理論預(yù)測值55.95s±0.07s較接近。

在1.0mol/L的NaCl溶液中，擴(kuò)散源個數(shù)n分別為10、20、30 和40 時，模擬流場的平均速度隨時間的變化關(guān)系如圖8所示。模擬流場的平均速度隨時間的變化先增大后減小，即相較于前期擴(kuò)散對流場速度的影響，Rayleigh 對流可有效加快流場速度；隨著對流結(jié)構(gòu)趨于穩(wěn)定，流場速度也開始下降。同時比較不同擴(kuò)散源個數(shù)的平均速度變化趨勢可知，隨著擴(kuò)散源個數(shù)的增加，速度的變化也會延遲。比較模擬得到的流場平均速度值與Thomas等[39]實驗中測得的流場速度1.82×10-4m/s 可知，模擬結(jié)果與實驗值具有相同的數(shù)量級。

2.3.2 鹽濃度對傳質(zhì)通量的影響

瞬時傳質(zhì)通量即單位時間通過單位截面積所傳遞物質(zhì)的質(zhì)量。由模擬結(jié)果得到的不同時刻CO2在流場中的濃度分布情況，可進(jìn)一步計算傳質(zhì)過程中不同時刻的瞬時傳質(zhì)通量Nins,t[kg/(m2·s)]，見式(22)。

圖9 給出了界面擴(kuò)散源數(shù)n=7 時CO2在純水和c0=1.0mol/L、2.0mol/L、3.0mol/L的咸水中瞬時傳質(zhì)通量Nins,t隨時間的變化關(guān)系。

觀察圖9 可知，對于不同鹽濃度的溶液，Nins,t隨時間增長的變化趨勢是一致的，都經(jīng)歷了先減小后增大再減小的變化。即前期氣液界面的傳質(zhì)以分子擴(kuò)散為主，由Fick定律可知，隨著CO2的不斷溶解，界面處的濃度梯度減小，傳質(zhì)速率變慢；而后隨著界面處與液相主體的密度差不斷增大，界面失穩(wěn)，高密度液體克服黏性力向下流動，并將液相主體的溶液擠壓到上方，促進(jìn)了液體表面更新，Rayleigh 對流結(jié)構(gòu)不斷向下和向兩側(cè)發(fā)展并相互誘導(dǎo)融合，使得對流結(jié)構(gòu)在溶液中快速擴(kuò)展，Nins,t也因此增大；最后，下邊界阻礙了Rayleigh對流結(jié)構(gòu)的發(fā)展，同時主體流場的濃度梯度差已經(jīng)較小，抑制了CO2在溶液中的擴(kuò)散，Nins,t開始減小。

圖6 n=7時不同鹽濃度和不同時刻溶質(zhì)的濃度分布

圖7 擴(kuò)散源個數(shù)n=1～40時Rayleigh對流開始時間的變化

圖8 n=10、20、30、40時流場平均速度隨時間變化關(guān)系

圖9 n=7時CO2在純水和不同鹽濃度咸水中瞬時傳質(zhì)通量Nins，t隨時間變化關(guān)系

將模擬結(jié)果與George 等[44]計算結(jié)果對比可知，兩者的瞬時傳質(zhì)通量的數(shù)量級［mg/(m2·s)］相同。同時，觀察圖9 中各個鹽濃度下Nins,t的最高點出現(xiàn)時間可知，純水和c0=1.0mol/L、2.0mol/L、3.0mol/L咸水中的出現(xiàn)時間分別為t=100s、125s、155s、200s左右，即鹽濃度的增加造成了Rayleigh對流發(fā)展的滯后，這與圖6中結(jié)果以及前人實驗得到的結(jié)論相吻合。同時觀察不同鹽濃度下Nins,t的變化曲線可知，鹽濃度的增加會降低CO2的瞬時傳質(zhì)通量。由此可知，鹽濃度的增加會抑制Rayleigh對流的發(fā)生和發(fā)展，使氣液傳質(zhì)速率總體降低。

3 結(jié)論

（1） 成功建立了用于咸水吸收CO2過程Rayleigh 對流模擬的LBM 模型，模擬得到的Rayleigh 對流結(jié)構(gòu)、對流發(fā)生時間、流場速度和瞬時傳質(zhì)通量的數(shù)量級與前人研究結(jié)果一致。

（2）與分子擴(kuò)散相比，Rayleigh 對流可有效強(qiáng)化咸水吸收CO2的氣液傳質(zhì)過程，且Rayleigh 對流結(jié)構(gòu)的發(fā)展規(guī)律與瞬時傳質(zhì)通量變化趨勢相關(guān)聯(lián)。

（3） 壁面阻力會抑制咸水吸收CO2過程Rayleigh對流的產(chǎn)生和發(fā)展。

（4）鹽濃度的增加會抑制Rayleigh對流的發(fā)生和發(fā)展，使氣液傳質(zhì)速率總體降低。

本文的模擬結(jié)果直觀展現(xiàn)了CO2封存于咸水層的過程中Rayleigh對流結(jié)構(gòu)對傳質(zhì)過程的促進(jìn)作用以及不同鹽濃度對傳質(zhì)效率的影響，對于實際的CO2封存過程有重要的指導(dǎo)意義。

符號說明

AI—— 截面面積，m2

Bi—— 速度系數(shù)，量綱為1

C—— 溶質(zhì)濃度，kg/m3

c—— 格子速度，量綱為1

c0—— 鹽濃度，mol/L

ci—— 離散速度，m/s

D—— 溶質(zhì)的擴(kuò)散系數(shù)，m2/s3

E—— 離子與CO2分子間靜電力作用勢能，J

EA—— CO2分子受到的總作用力勢能，J

e—— 單位電荷的帶電量，1.602177×10-19C

Fp—— 單位體積的體積力，N/m3

Fw—— 單位體積的壁面阻力，N/m3

f，g—— 粒子分布函數(shù)

g—— 重力加速度，m/s2

H—— 計算區(qū)域高度，mm

L—— 計算區(qū)域長度，mm

Nins,t—— 瞬時傳質(zhì)通量，kg/(m2·s）

n—— 溶質(zhì)擴(kuò)散源個數(shù)

tonset—— 對流開始時間，s

Δt—— 單位時間步長，s

u—— 宏觀速度，m/s

V—— 溶液體積，m3

Δx—— 單位網(wǎng)格步長，m

α—— CO2極化率，2.911×10-24C·m2/V

δ—— 壁面位置，mm

ε0—— 真空介電常數(shù)，8.854187817×10-12F/m

μ—— 溶液動力黏度，Pa·s

ν—— 液體運(yùn)動黏度，2.911×10-24m2/s

ρ—— 液體宏觀密度，kg/m3

Δρ—— 液體密度差，kg/m3

τ—— 松弛時間，量綱為1

ω—— 權(quán)重系數(shù)，量綱為1

下角標(biāo)

A—— 溶劑

B—— 溶質(zhì)

i—— 粒子運(yùn)動方向，i=1,…,9

ins—— 瞬時