多孔介質通道內非混相驅替過程的格子Boltzmann方法模擬

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

專題報道

多孔介質通道內非混相驅替過程的格子Boltzmann方法模擬

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

（華東理工大學 化工學院 化學工程聯合國家重點實驗室，上海 200237）

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

Shan-Chen偽勢格子Boltzmann方法；驅替；二氧化碳

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

格子Boltzmann方法（LBM）具有清晰的粒子運動圖像，在滲流等方面得到了廣泛應用。Kang等［6］利用LBM模型模擬了兩相不互溶流體在重力作用下沿固壁鋪展動態過程；Li等［7］利用LBM偽勢模型，通過改進液固作用力，實現了較大密度比兩組分流體對壁面潤濕的模擬。

本工作基于多組分LBM模型、耦合相間作用力，研究了雷諾數（Re）、毛細數（Ca）和壁面潤濕性對非混相驅替流動中黏性指進的影響，探索了CO2泡沫在多孔介質內流動的物理本質，對CO2驅油法提高油氣采收率（CO2-EOR）技術的推廣應用有著極其重要的現實意義。

1 LBM理論

LBGK模型是應用最為廣泛的LBE（Lattice Boltzmann equation）模型，本工作以D2Q9 LBGK模型對多孔介質通道內的滲流過程進行數值模擬。在LBGK模型中，粒子的運動表現為碰撞和遷移兩個步驟，如式（1）所示：

方程左邊代表遷移過程，右邊為碰撞過程。其中，fασ（x，t）是第 σ 種組分在 x空間位置、t時刻的 α 速度方向上的密度分布函數，fασ（eq）（x，t）是其對應的平衡態密度分布函數，τσ是σ組分的無量綱松弛時間，與該組分的動力學黏度有關（vσ=（τσ- 0.5）δt/3）。

圖1給出了D2Q9速度離散模型示意圖，其各個方向的離散速度見式（2）：

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

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

2 模擬設置

2.1 計算區域及邊界條件設置

圖2a為驅替流動示意圖和黏性指進考察參數（D/L）的定義。孔道的長與寬分別為lx，ly，計算區域網格數為200×50。初始狀態，孔道內充滿被驅替流體，始終處于靜止狀態；驅替開始，驅替流體以速度u0進入孔道左側，右側為流體出口。上、下側為通道壁面。在地藏條件下，CO2泡沫實質是超臨界CO2水乳液［10-11］，相對密度在0.40～0.85之間，原油相對密度在0.75～1.00之間［12-13］，兩相密度比在0.88～3.00范圍內，本工作模擬中取兩相流體密度比為1.00。

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

圖2 驅替流動示意圖（a）和黏性指進考察參數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定律給出了靜止液滴的內、外壓力差與靜止液滴半徑之間的關系，見式（4）。

pin-pout= γ/r （4）

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

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

由此，可通過改變液-液相互作用力參數Gx值的大小來模擬界面張力對兩相流動的影響。

2.3 接觸角演化

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

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

圖5給出了三種典型潤濕情況下，液滴在固壁上的鋪展情況。從圖5可見，當Gw1= 0.02時，θ=120.51°，為非潤濕流體，液滴在固壁上呈收縮聚集；當Gw1= 0時，θ= 90.00°，液滴呈半球狀，固壁對液滴無作用力；當Gw1= -0.02時，θ= 45.34°，為潤濕性流體，液滴鋪展在固壁上。

圖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.

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

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

由圖6可見，θ隨Gw1的增大線性增加，線性擬合關聯式見式（7）。

根據結果，可通過調節Gw1值模擬對固壁不同潤濕性的情況下的驅替流動。

3 孔道驅流黏性指進的影響因素

根據界面形狀，定義了無量綱量界面結構參數D/L度量流動中黏性指進的程度，D/L值越大，兩相界面越彎曲，表明黏性指進情況越嚴重。同時，定義驅替流體流動界面抵達出口處時驅替流體所占孔道總面積比率為驅替波及面積效率（Se），用于度量驅替效果，Se越大，波及面積越大，驅替效果越好。

3.1 雷諾數的影響

為考察Re對孔道內兩相驅替流動發生的黏性指進，考察了不同Re下的驅流情況。模擬中不考慮流-固間的相互作用，假設Gw1=Gw2= 0。驅替流體的注入速度在0.01～0.05之間變化，此時對應Re= 3.333～16.667。Gx等其他模型參數不變。

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

圖7 不同Re下兩相界面隨時間的演化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時，不同Re下驅替流動相界面的D/L值時發現（圖8），D/L隨Re的增加而增大。當Re從3.333增加至16.667時，Se值減小了15.80%（圖8b）。由此可知，為保證獲得較大的驅替效率，必須根據物性提供合理的注入速度，否則黏性指進現象不可避免。

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

3.2 毛細數的影響

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

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

圖9 不同Ca下相界面隨時間的演化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對D/L（a）和Se（b）的影響Fig.10 Effect of Ca on D/L (a) and Se(b).

3.3 潤濕性的影響

圖11為不同液固接觸角條件下驅替時間為t=3 000時的模擬結果。由圖11可知，驅替流體相對固體壁面潤濕程度越佳，則相界面越平整，黏性指進程度越弱；若兩相間非潤濕性越強，則相界面越彎曲，黏性指進程度越嚴重。

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

圖11 不同壁面潤濕性條件下的驅替情況Fig.11 Displacement simulations for different surface wet ability.

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

4 結論

1）利用Shan-Chen偽勢格子Boltzmann模型模擬了超臨界CO2泡沫驅油流動，當Re、Ca較大時，均會導致黏性指進現象顯著發展，降低驅替效率。

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

3）非潤濕流體極易導致黏性指進現象，且與非潤濕程度強相關，當接觸角從57.77°增加到111.73°時，Se減小了18.35%。在進行兩相驅替流動時，須綜合物性及驅替流體與固壁的潤濕情況，確定驅替流體的注入條件，以避免黏性指進現象，提高驅替效率。

符 號 說 明

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

a 液滴的高度

b 平板上的鋪展寬度

Ca 毛細數

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

Error相鄰兩次全局速度相對誤差

eα粒子離散速度

Fσ，Ffσ，Fwσ，Feσ各作用力，液-液相互作用力，液-固

相互作用力，外部作用力

fασ，fασ(eq)分布函數，平衡態分布函數

Gx，Gw，Gw1，Gσσ液-液相互作用參數，液-固相互作用參數，組分1對固壁作用力參數，各組分間相互作用力

lx，ly計算域尺寸

p，pin，pout壓力，液滴內壓力，液滴外壓力

Re 雷諾數

r 液滴半徑

Se 驅替波及面積效率

t 流動時間

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

態速度，混合流體速度

x，y 驅替方向及流動法向方向

γ 界面張力

δt，δx時間步長，空間步長

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

υ 黏度系數

ρσ密度

σ，σ 第σ及σ種組分

τ 松弛時間

ψσ（x） 有效密度

ωα權重系數

［1］ 劉朝全，姜學峰. 2016年國內外油氣行業發展報告［M］.北京：石油工業出版社，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］ 柳瑤閣. 多孔介質重構與指進模擬研究［D］.合肥：中國科學技術大學，2014.

［5］ 陳挺. 稠油聚/表復合驅波及系數和驅油效率影響因素及調控方法研究［D］.青島：中國石油大學（華東），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］ 秦積舜，韓海水，劉曉蕾. 美國 CO2驅油技術應用及啟示［J］.石油勘探與開發，2015，42（2）：209-216.

［12］ 蘇玉亮. 油藏驅替機理［M］.北京：石油工業出版社，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］ 李盂濤，單文文，劉先貴. 超臨界二氧化碳混相驅油機理實驗研究［J］.石油學報，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—），男，河南省沁陽市人，碩士生，電話 13127581236，電郵 liujwecust@163.com。聯系人：宗原，電話 021-64253042，電郵 zongyuan@ecust.edu.cn。