劉佳威，許志美，宗 原，劉 濤

專題報(bào)道

多孔介質(zhì)通道內(nèi)非混相驅(qū)替過(guò)程的格子Boltzmann方法模擬

劉佳威，許志美，宗 原，劉 濤

（華東理工大學(xué) 化工學(xué)院 化學(xué)工程聯(lián)合國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200237）

針對(duì)CO2泡沫相與原油體系，基于Shan-Chen偽勢(shì)格子Boltzmann模型，耦合液-液、液-固間相互作用力與界面張力、液固潤(rùn)濕性之間的關(guān)系，對(duì)多孔介質(zhì)通道內(nèi)CO2泡沫驅(qū)流動(dòng)規(guī)律進(jìn)行數(shù)值模擬，研究了雷諾數(shù)（Re）、毛細(xì)數(shù)（Ca）和壁面潤(rùn)濕性對(duì)非混相兩相驅(qū)替流動(dòng)中黏性指進(jìn)的影響。實(shí)驗(yàn)結(jié)果表明，隨Re和Ca的增大，黏性指進(jìn)顯著發(fā)展，驅(qū)替效率降低；壁面潤(rùn)濕時(shí)，無(wú)明顯指進(jìn)現(xiàn)象，而當(dāng)壁面非潤(rùn)濕時(shí)，黏性指進(jìn)較為明顯，且隨非潤(rùn)濕性程度的增加而愈發(fā)嚴(yán)重。在進(jìn)行兩相驅(qū)替流動(dòng)時(shí)，需綜合物性及驅(qū)替流體與固壁的潤(rùn)濕情況，確定驅(qū)替流體的注入條件，以提高驅(qū)替效率。

Shan-Chen偽勢(shì)格子Boltzmann方法；驅(qū)替；二氧化碳

我國(guó)原油的需求量以每年2%～3%的速度增長(zhǎng)［1］，低滲透油氣資源的勘探開(kāi)發(fā)被提上日程。其中，CO2泡沫驅(qū)替技術(shù)成為提高低滲透油藏的重要技術(shù)之一［2-3］。但在驅(qū)替過(guò)程中極易產(chǎn)生黏性指進(jìn)，影響驅(qū)油效率。研究發(fā)現(xiàn)黏性指進(jìn)現(xiàn)象與兩相黏度比、驅(qū)替注入速度、壁面潤(rùn)濕性和原油與驅(qū)替相的界面張力等因素［4-5］相關(guān)。

格子Boltzmann方法（LBM）具有清晰的粒子運(yùn)動(dòng)圖像，在滲流等方面得到了廣泛應(yīng)用。Kang等［6］利用LBM模型模擬了兩相不互溶流體在重力作用下沿固壁鋪展動(dòng)態(tài)過(guò)程；Li等［7］利用LBM偽勢(shì)模型，通過(guò)改進(jìn)液固作用力，實(shí)現(xiàn)了較大密度比兩組分流體對(duì)壁面潤(rùn)濕的模擬。

本工作基于多組分LBM模型、耦合相間作用力，研究了雷諾數(shù)（Re）、毛細(xì)數(shù)（Ca）和壁面潤(rùn)濕性對(duì)非混相驅(qū)替流動(dòng)中黏性指進(jìn)的影響，探索了CO2泡沫在多孔介質(zhì)內(nèi)流動(dòng)的物理本質(zhì)，對(duì)CO2驅(qū)油法提高油氣采收率（CO2-EOR）技術(shù)的推廣應(yīng)用有著極其重要的現(xiàn)實(shí)意義。

1 LBM理論

LBGK模型是應(yīng)用最為廣泛的LBE（Lattice Boltzmann equation）模型，本工作以D2Q9 LBGK模型對(duì)多孔介質(zhì)通道內(nèi)的滲流過(guò)程進(jìn)行數(shù)值模擬。在LBGK模型中，粒子的運(yùn)動(dòng)表現(xiàn)為碰撞和遷移兩個(gè)步驟，如式（1）所示：

方程左邊代表遷移過(guò)程，右邊為碰撞過(guò)程。其中，fασ（x，t）是第 σ 種組分在 x空間位置、t時(shí)刻的 α 速度方向上的密度分布函數(shù)，fασ（eq）（x，t）是其對(duì)應(yīng)的平衡態(tài)密度分布函數(shù)，τσ是σ組分的無(wú)量綱松弛時(shí)間，與該組分的動(dòng)力學(xué)黏度有關(guān)（vσ=（τσ- 0.5）δt/3）。

圖1給出了D2Q9速度離散模型示意圖，其各個(gè)方向的離散速度見(jiàn)式（2）：

圖1 D2Q9模型離散速度設(shè)定Fig.1 Discrete velocity set of D2Q9 model.ei：the discrete-velocity of each direction，i = 0 - 8.

此外，在D2Q9模型中，流體混合物的總壓力見(jiàn)式（3）［9］。

2 模擬設(shè)置

2.1 計(jì)算區(qū)域及邊界條件設(shè)置

圖2a為驅(qū)替流動(dòng)示意圖和黏性指進(jìn)考察參數(shù)（D/L）的定義?？椎赖拈L(zhǎng)與寬分別為lx，ly，計(jì)算區(qū)域網(wǎng)格數(shù)為200×50。初始狀態(tài)，孔道內(nèi)充滿被驅(qū)替流體，始終處于靜止?fàn)顟B(tài)；驅(qū)替開(kāi)始，驅(qū)替流體以速度u0進(jìn)入孔道左側(cè)，右側(cè)為流體出口。上、下側(cè)為通道壁面。在地藏條件下，CO2泡沫實(shí)質(zhì)是超臨界CO2水乳液［10-11］，相對(duì)密度在0.40～0.85之間，原油相對(duì)密度在0.75～1.00之間［12-13］，兩相密度比在0.88～3.00范圍內(nèi)，本工作模擬中取兩相流體密度比為1.00。

對(duì)上、下固體壁面，采用反彈格式處理實(shí)現(xiàn)非滑移邊界條件，進(jìn)口邊界采用Zou-He速度邊界條件［14］處理，出口邊界為充分發(fā)展邊界條件。計(jì)算初始，驅(qū)替相與被驅(qū)替的密度分布函數(shù)均為以初速度、初始態(tài)密度計(jì)算的fσeq（x，t），隨后，進(jìn)行遷移、碰撞步驟迭代計(jì)算直至結(jié)束。

圖2 驅(qū)替流動(dòng)示意圖（a）和黏性指進(jìn)考察參數(shù)D/L的定義（b）Fig.2 Schematic diagram of displacement flow(a) and definition of viscous fingering parameter D/L(b).x，y：the direction along and normal to the displacement path，respectively；lx，ly：two dimensions of the computational domain；D：the extension length of the finger；L：the half width of the channel.

2.2 界面張力演化

Laplace-Yang定律給出了靜止液滴的內(nèi)、外壓力差與靜止液滴半徑之間的關(guān)系，見(jiàn)式（4）。

pin-pout= γ/r （4）

圖3a為不同Gx下 （pin- pout）～ 1/r的關(guān)系曲線。

從圖3a可看出，液滴內(nèi)、外壓差與液滴半徑倒數(shù)呈線性關(guān)系，與Laplace-Yang定律相符。根據(jù)Laplace-Yang定律，整理圖中數(shù)據(jù)可得兩相界面張力γ與Gx間線性關(guān)系，如圖3b所示。經(jīng)擬合可得式（5）：

由此，可通過(guò)改變液-液相互作用力參數(shù)Gx值的大小來(lái)模擬界面張力對(duì)兩相流動(dòng)的影響。

2.3 接觸角演化

接觸角θ用于表征流體對(duì)固體壁面的親、疏水性質(zhì)，θ＜90°，θ＞90°分別對(duì)應(yīng)于潤(rùn)濕性、疏水性流體。為了建立流-固作用力參數(shù)Gw與接觸角的關(guān)系，模擬組分1在固體壁面上的潤(rùn)濕鋪展過(guò)程。計(jì)算域網(wǎng)格設(shè)置為lxly= 100×100，初始設(shè)定底部中間靜置一個(gè)半徑為r的圓形液滴（組分1），其余部分被流體2占據(jù)。密度設(shè)定ρ1=ρ2= 2.0，兩相對(duì)固壁作用力參數(shù)Gw1= -Gw2，液滴在液-液作用力與流-固作用力的共同作用下自由演化（見(jiàn)圖4）。液滴穩(wěn)定后，接觸角θ可由液滴高度a與底部寬度b［15］獲得，見(jiàn)式（6）。

圖4 接觸角計(jì)算示意圖Fig.4 Schematic diagram of contact angle computation.a：the height of the drop；b：the width of the drop；θ：the contact angle.

圖5給出了三種典型潤(rùn)濕情況下，液滴在固壁上的鋪展情況。從圖5可見(jiàn)，當(dāng)Gw1= 0.02時(shí)，θ=120.51°，為非潤(rùn)濕流體，液滴在固壁上呈收縮聚集；當(dāng)Gw1= 0時(shí)，θ= 90.00°，液滴呈半球狀，固壁對(duì)液滴無(wú)作用力；當(dāng)Gw1= -0.02時(shí)，θ= 45.34°，為潤(rùn)濕性流體，液滴鋪展在固壁上。

圖5 不互溶兩相流體與壁面接觸角的模擬Fig.5 Simulation of different contact angles for two immiscible fluids.(a) Gw1 = 0.02，θ = 120.51°，non-wetting surface；(b) Gw1 = 0，θ = 90.00°，mediate surface；(c) Gw1 = -0.02，θ = 45.34°，wetting surface.Gw1：Liquid-solid interaction parameter.

基于此，可獲得接觸角θ與流-固作用力參數(shù)Gw1之間的關(guān)系曲線（圖6）。

圖6 液-固作用力參數(shù)與接觸角的關(guān)系曲線Fig.6 Relationships between contact angle and liquid-solid interaction parameter.

由圖6可見(jiàn)，θ隨Gw1的增大線性增加，線性擬合關(guān)聯(lián)式見(jiàn)式（7）。

根據(jù)結(jié)果，可通過(guò)調(diào)節(jié)Gw1值模擬對(duì)固壁不同潤(rùn)濕性的情況下的驅(qū)替流動(dòng)。

3 孔道驅(qū)流黏性指進(jìn)的影響因素

根據(jù)界面形狀，定義了無(wú)量綱量界面結(jié)構(gòu)參數(shù)D/L度量流動(dòng)中黏性指進(jìn)的程度，D/L值越大，兩相界面越彎曲，表明黏性指進(jìn)情況越嚴(yán)重。同時(shí)，定義驅(qū)替流體流動(dòng)界面抵達(dá)出口處時(shí)驅(qū)替流體所占孔道總面積比率為驅(qū)替波及面積效率（Se），用于度量驅(qū)替效果，Se越大，波及面積越大，驅(qū)替效果越好。

3.1 雷諾數(shù)的影響

為考察Re對(duì)孔道內(nèi)兩相驅(qū)替流動(dòng)發(fā)生的黏性指進(jìn)，考察了不同Re下的驅(qū)流情況。模擬中不考慮流-固間的相互作用，假設(shè)Gw1=Gw2= 0。驅(qū)替流體的注入速度在0.01～0.05之間變化，此時(shí)對(duì)應(yīng)Re= 3.333～16.667。Gx等其他模型參數(shù)不變。

圖7 為不同Re下兩相界面隨時(shí)間的演化情況。圖7中實(shí)線所示為流體流動(dòng)過(guò)程中經(jīng)Δt時(shí)間間隔后兩相界面輪廓。從圖7可見(jiàn)，低Re下，兩相驅(qū)替界面較平，基本沒(méi)有指進(jìn)現(xiàn)象形成（如圖7a和b）。當(dāng)Re增加至10.00（7c，d，e），兩相界面不再平坦，“手指”逐漸突顯，黏性指進(jìn)現(xiàn)象明顯。

圖7 不同Re下兩相界面隨時(shí)間的演化Fig.7 Evolution of the interfaces with time for different Re.(a) Re=3.333，Δt = 2 000；(b) Re = 6.667，Δt = 1 000；(c) Re = 10.00，Δt = 600；(d) Re = 13.333，Δt = 600；(e) Re = 16.667，Δt = 400 Δt：time interval.

分析t= 3 000時(shí)，不同Re下驅(qū)替流動(dòng)相界面的D/L值時(shí)發(fā)現(xiàn)（圖8），D/L隨Re的增加而增大。當(dāng)Re從3.333增加至16.667時(shí)，Se值減小了15.80%（圖8b）。由此可知，為保證獲得較大的驅(qū)替效率，必須根據(jù)物性提供合理的注入速度，否則黏性指進(jìn)現(xiàn)象不可避免。

圖8 Re對(duì)D/L（a）和驅(qū)替波及面積效率（Se）（b）的影響Fig.8 Effection of Re on D/L(a) and sweeping efficiency(Se)(b).

3.2 毛細(xì)數(shù)的影響

Ca表示黏性力相對(duì)表面張力的大小，Ca=υρu/γ。根據(jù)式（5），選取不同Gx值考察Ca對(duì)驅(qū)流的影響。驅(qū)替流體速度設(shè)為0.02，Gw保持不變。圖9分別列出了不同Ca下孔道中兩相驅(qū)替界面隨時(shí)間的演化過(guò)程。圖中相鄰界面經(jīng)過(guò)時(shí)間間隔Δt為1 000。從圖9可見(jiàn)，在低Ca情況下，指進(jìn)情況變化不大。Ca增至0.030 62（即γ =0.196 0）時(shí)，才可明顯觀察到指進(jìn)現(xiàn)象。且隨驅(qū)替流動(dòng)時(shí)間的增加，“手指”現(xiàn)象逐漸凸顯，顯示出相界面與壁面接觸點(diǎn)的滯后現(xiàn)象。

圖10給出了t= 3 000時(shí)不同Ca下驅(qū)替流動(dòng)相界面處的D/L值。從圖10 可知，D/L值隨Ca的增加而提高，且增加趨勢(shì)逐漸加劇。Se隨Ca的增大變化不大，Ca從0.015 16增加到0.030 62，Se僅減小約6.07%?；煜囹?qū)替與非混相驅(qū)替過(guò)程的作用機(jī)理不同?；煜囹?qū)替時(shí)，低界面張力有助于抑制指進(jìn)現(xiàn)象，提高驅(qū)替效率［16］，而對(duì)于非混相驅(qū)替而言，兩相間作用力增加，Ca減小，則有利于提高波及效率，這與文獻(xiàn)［17］中的結(jié)論一致。

圖9 不同Ca下相界面隨時(shí)間的演化Fig.9 Time evolution of the interfaces for different Capillary number(Ca).Ca：(a) 0.015 16；(b) 0.016 84；(c) 0.019 34；(d) 0.023 47；(e) 0.030 62

圖10 Ca對(duì)D/L（a）和Se（b）的影響Fig.10 Effect of Ca on D/L (a) and Se(b).

3.3 潤(rùn)濕性的影響

圖11為不同液固接觸角條件下驅(qū)替時(shí)間為t=3 000時(shí)的模擬結(jié)果。由圖11可知，驅(qū)替流體相對(duì)固體壁面潤(rùn)濕程度越佳，則相界面越平整，黏性指進(jìn)程度越弱；若兩相間非潤(rùn)濕性越強(qiáng)，則相界面越彎曲，黏性指進(jìn)程度越嚴(yán)重。

對(duì)圖11中的驅(qū)替現(xiàn)象進(jìn)行整理，可得到D/L值和隨驅(qū)替流體與固體壁面接觸角θ的變化關(guān)系（圖12a）。從圖12可看出，在50°lt;θlt; 80°范圍內(nèi)，即對(duì)潤(rùn)濕液體而言，黏性指進(jìn)均較弱，且黏性指進(jìn)基本不受壁面潤(rùn)濕性影響。隨壁面非潤(rùn)濕性增強(qiáng)，黏性指進(jìn)現(xiàn)象加劇。圖12b為Se與固壁潤(rùn)濕條件的關(guān)系。從圖12b可看出，潤(rùn)濕較好時(shí)（50°lt;θlt;80°）驅(qū)替效率較高，隨非潤(rùn)濕性的增加，Se急劇減小，且在θ=110°處趨于穩(wěn)定。當(dāng)θ從57.77°增加到111.73°時(shí)，Se減小了18.35%。

圖11 不同壁面潤(rùn)濕性條件下的驅(qū)替情況Fig.11 Displacement simulations for different surface wet ability.

圖12 接觸角對(duì)D/L（a）和Se （b）的影響Fig.12 Effection of contact angle on D/L(a) and Se(b).

4 結(jié)論

1）利用Shan-Chen偽勢(shì)格子Boltzmann模型模擬了超臨界CO2泡沫驅(qū)油流動(dòng)，當(dāng)Re、Ca較大時(shí)，均會(huì)導(dǎo)致黏性指進(jìn)現(xiàn)象顯著發(fā)展，降低驅(qū)替效率。

2）隨Re從3.333增加到16.667，Se減小了15.80%；Ca從0.015 16增加到0.030 62，Se減小了6.07%。

3）非潤(rùn)濕流體極易導(dǎo)致黏性指進(jìn)現(xiàn)象，且與非潤(rùn)濕程度強(qiáng)相關(guān)，當(dāng)接觸角從57.77°增加到111.73°時(shí)，Se減小了18.35%。在進(jìn)行兩相驅(qū)替流動(dòng)時(shí)，須綜合物性及驅(qū)替流體與固壁的潤(rùn)濕情況，確定驅(qū)替流體的注入條件，以避免黏性指進(jìn)現(xiàn)象，提高驅(qū)替效率。

符 號(hào) 說(shuō) 明

（除接觸角外其他參數(shù)單位均基于格子單位）

a 液滴的高度

b 平板上的鋪展寬度

Ca 毛細(xì)數(shù)

c，cs格子速度，格子聲速

Error相鄰兩次全局速度相對(duì)誤差

eα粒子離散速度

Fσ，F(xiàn)fσ，F(xiàn)wσ，F(xiàn)eσ各作用力，液-液相互作用力，液-固

相互作用力，外部作用力

fασ，fασ(eq)分布函數(shù)，平衡態(tài)分布函數(shù)

Gx，Gw，Gw1，Gσσ液-液相互作用參數(shù)，液-固相互作用參數(shù)，組分1對(duì)固壁作用力參數(shù)，各組分間相互作用力

lx，ly計(jì)算域尺寸

p，pin，pout壓力，液滴內(nèi)壓力，液滴外壓力

Re 雷諾數(shù)

r 液滴半徑

Se 驅(qū)替波及面積效率

t 流動(dòng)時(shí)間

u，u0，uσ，uσeq，u′ 速度，初速度，各組分速度，平衡

態(tài)速度，混合流體速度

x，y 驅(qū)替方向及流動(dòng)法向方向

γ 界面張力

δt，δx時(shí)間步長(zhǎng)，空間步長(zhǎng)

ε 相鄰兩次全局速度相對(duì)誤差θ 接觸角，°

υ 黏度系數(shù)

ρσ密度

σ，σ 第σ及σ種組分

τ 松弛時(shí)間

ψσ（x） 有效密度

ωα權(quán)重系數(shù)

［1］ 劉朝全，姜學(xué)峰. 2016年國(guó)內(nèi)外油氣行業(yè)發(fā)展報(bào)告［M］.北京：石油工業(yè)出版社，2017.

［2］ Farajzadeh R，Andrianov A，Bruining H，et al. Comparative study of CO2and N2foams in porous media at low and high pressure-temperatures［J］.Ind Eng Chem Res，2009，48（9）：4542-4552.

［3］ Emadi A，Jamiolahmady M，Sohrabi M，et al. Visualization of oil recovery by CO2-foam injection；effect of oil viscosity and gas type［C］//SPE Improved Oil Recovery Symposium. Society of Petroleum Engineers，2012.

［4］ 柳瑤閣. 多孔介質(zhì)重構(gòu)與指進(jìn)模擬研究［D］.合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2014.

［5］ 陳挺. 稠油聚/表復(fù)合驅(qū)波及系數(shù)和驅(qū)油效率影響因素及調(diào)控方法研究［D］.青島：中國(guó)石油大學(xué)（華東），2014.

［6］ Kang Qinjun，Zhang Dongxiao，Chen Shiyi. Displacement of a two-dimensional immiscible droplet in a channel［J］.Phys Fluids，2002，14（9）：3203-3214.

［7］ Li Qing，Luo Kun，Kang Qinjun，et al. Contact angles in the pseudopotential lattice Boltzmann modeling of wetting［J］.Phys Rev E，2014，90（5）：053301.

［8］ Shan X，Doolen G. Multicomponent lattice-Boltzmann model with interparticle interaction［J］.J Stat Phys，1995，81（1/2）：379-393.

［9］ Shan X，Doolen G. Diffusion in a multicomponent lattice Boltzmann equation model［J］.Phys Rev，1996，54（4）：3614-3620.

［10］ Saghir M Z，Chaalal O，Islam M R. Numerical and experimental modeling of viscous fingering during liquid-liquid miscible displacement［J］.J Pet Sci Eng，2000，26（1）：253-262.

［11］ 秦積舜，韓海水，劉曉蕾. 美國(guó) CO2驅(qū)油技術(shù)應(yīng)用及啟示［J］.石油勘探與開(kāi)發(fā)，2015，42（2）：209-216.

［12］ 蘇玉亮. 油藏驅(qū)替機(jī)理［M］.北京：石油工業(yè)出版社，2009.

［13］ Anwar S，Carroll J J. Carbon dioxide thermodynamic properties handbook：Covering temperatures from -20 degrees to 250 degrees celcius and pressures up to 1 000 bar［M］.Alberta：John Wiley amp; Sons，2011：15-50.

［14］ Zou Qisu，He Xiaoyi. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model［J］.Phys Fluids，1997，9（6）：1591-1598.

［15］ Dullien F A L. Porous media：Fluid transport and pore structure［M］.Ontario：Academic Press，2012.

［16］ 李盂濤，單文文，劉先貴. 超臨界二氧化碳混相驅(qū)油機(jī)理實(shí)驗(yàn)研究［J］.石油學(xué)報(bào)，2006，27（3）：80-83.

［17］ Chin J. Lattice Boltzmann simulation of the flow of binary immiscible fluids with different viscosities using the Shan-Chen microscopic interaction model［J］.Philos Trans，2002，360（1792）：547-558.

（編輯 平春霞）

Lattice Boltzmann simulation for CO2foam flow in porous media

Liu Jiawei，Xu Zhimei，Zong Yuan，Liu Tao
（Chemical Engineering Department，State Key Laboratory of Chemical Engineering，East China University of Science and Technology，Shanghai 200237，China）

The relationship between liquid-liquid interaction，solid-liquid interaction and surface tension，surface wettability is coupled with Shan-Chen pseudopotential lattice Boltzmann model，and the flow phenomena of carbon dioxide foam through porous media is simulated. The effects of Reynolds number(Re)，capillary number(Ca) and surface wettability on viscous fingering is investigated.The simulation has demonstrated that the viscous fingering phenomena will be enhanced significantly with the increase ofReandCa. While，no obvious viscous fingering can be observed when the surface is wetting to the displacing fluid. On the contrary，the viscous fingering phenomena is intensified when the surface is non-wetting to the displacing fluid. It’s illustrated that the injection conditions should be set according to the fluids properties and wettability to ensure high displacement efficiency.

Shen-Chen pseudopotential lattice Boltzmann method；displacement；carbon dioxide

1000-8144（2017）11-1347-08

TQ 021.1

A

10.3969/j.issn.1000-8144.2017.11.001

2017-03-10；［修改稿日期］2017-08-07。

劉佳威（1991—），男，河南省沁陽(yáng)市人，碩士生，電話 13127581236，電郵 liujwecust@163.com。聯(lián)系人：宗原，電話 021-64253042，電郵 zongyuan@ecust.edu.cn。