多孔介质通道内非混相驱替过程的格子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。