STUDY ON SNAP-OFF MECHANISM AND SIMULATION DURING GAS-LIQUID IMMISCIBLE DISPLACEMENT
-
摘要: 研究气液非混相驱替过程中的相界面卡断机理及其影响因素在气驱, 气水交替及泡沫驱等提高油气采收率领域具有重要意义. 本文在原始伪势格子玻尔兹曼模型的基础上, 改进流体-流体作用力格式, 添加流-固作用力, 耦合RK状态方程, 并采用精确差分方法将外力添加到LBM框架中. 通过校准模型的热力学一致性以及模拟测试界面张力, 静态平衡接触角及液相在角隅的滞留等一系列两相体系验证模型的准确性. 基于改进的伪势格子玻尔兹曼模型, 在孔-喉-孔系统中开展气液非混相驱替模拟, 结果表明: 卡断现象与驱替压差, 孔喉长度比及孔喉宽度比有关, 只有当驱替压差处于一定范围内时, 气液两相驱替过程中才会发生卡断现象; 当驱替压差大于临界驱替压差上限时, 即使达到了经典静态准则所预测的卡断条件, 卡断也会被抑制; 当驱替压差小于临界驱替压差下限时, 无法克服毛管“钉扎”作用, 形成无效驱替. 对于固定孔喉宽度比的孔-喉-孔结构, 随着孔喉长度比的增大, 发生卡断现象的驱替压差范围增大; 对于固定孔喉长度比的孔-喉-孔结构, 随着孔喉宽度比的减小, 发生卡断现象的驱替压差范围增大.Abstract: Studying the mechanism of phase interface snap-off during gas liquid immiscible displacement and its influencing factors have great significance in the field of enhanced oil and gas recovery such as gas driving, gas water alternation and foam driving. In this work, based on the original pseudopotential lattice Boltzmann model, we improved the fluid-fluid force scheme, added the fluid-solid force, coupled the Redlich-Kwong (RK) equation of state, and used the exact difference method (EDM) to add the external forces to the LBM framework. As well as verified the accuracy of the model by calibrating the thermodynamic consistency of the model and simulating a series of two phase systems such as testing the interfacial tension, static equilibrium contact angle and retention of the liquid phase at the corner. Based on the modified pseudopotential lattice Boltzmann model, we have carried out gas-liquid immiscible displacement simulations in a pore-throat-pore system, and the results have shown that: the snap-off phenomenon is related to the displacement pressure difference, pore-throat length ratio and pore-throat width ratio, and the snap-off phenomenon occurs only when the displacement pressure difference is within a certain range. When the displacement pressure difference is larger than the upper limit of the critical displacement pressure difference, the snap-off will be inhibited even if the snap-off condition predicted by the classical static rule has been reached; When the displacement pressure difference is less than the lower limit of the critical displacement pressure difference, it cannot overcome the "pinning" effect of the capillary tube and results in ineffective displacement. For the pore-throat structure with constant pore-throat width ratio, the displacement pressure difference range in which the snap-off phenomenon occurs increases as the pore-throat length ratio increases; For the pore-throat structure with constant pore-throat length ratio, the displacement pressure difference range in which the snap-off phenomenon occurs increases as the pore-throat width ratio decreases.
-
引 言
复杂多孔介质中的非混相驱替在两相流动中占据重要地位, 并在气驱, 气水交替及泡沫驱等提高油气采收率领域具有重要研究意义[1-5]. 在气驱提高油气采收率过程中, 连续气相注入储层后, 将会形成复杂的气液非混相驱替模式; 注入储层中的气相一般为非润湿相, 在非润湿相驱替润湿相的过程中, 孔喉壁面上将会滞留一层液膜, 液膜的存在可能会造成驱替过程中发生卡断甚至水锁现象. 因此, 在孔隙尺度研究气液非混相驱替过程中的卡断机理有利于人们认识微观孔喉中的气液流动状态.
研究人员通过理论研究, 实验及数值模拟等方法, 对气液两相驱替过程中的卡断机理展开了大量研究工作. 在理论方面, Roof[6]基于孔喉毛管压力不平衡这一机制, 建立了预测卡断模型的静态准则, 认为当最小喉道半径小于孔隙半径的一半时, 润湿相将沿着孔隙壁面发生回流, 也即导致非润湿相发生卡断. Gauglitz等[7]给出了收缩圆柱形孔喉系统的卡断准则, 认为毛管数Ca对卡断发生的时间有一定影响, 并且对于不同宽度比的孔喉系统影响规律不同. Ransohoff等[8]给出了方形喉道截面的以及三角形喉道截面的静态卡断准则, 并认为在非圆形截面的孔喉中更容易发生卡断现象. 然而, 上述卡断准则并没有考虑动态参数对卡断影响, Tsai和Miksis[9]通过数值方法发现毛管数过大或者过小时, 气泡不会发生卡断现象. Deng等[10-11]扩展了圆形及非圆形缩颈孔喉系统的静态卡断准则, 发现之前的静态准则对卡断的预测过于保守, 同时发现即使喉道系统达到了卡断发生的静态准则, 高毛管数也会抑制卡断的发生. 在实验方面, 随着微流控模型的发展, 大量学者采用微流控模型对两相流动的卡断现象进行可视化实验, Tian等[12]利用2D微模型研究了气水流动过程中的卡断现象, 并验证了Roof等提出的静态准则. Cha等[13]利用2D矩形截面微模型研究了气泡卡断机理, 并与提出的卡断理论模型得到了较好的一致性. Tetteh等[14]利用微流控模型观察到了低矿化度水驱过程中的油相卡断. Wu等[15]利用实验研究非缩颈孔喉的气泡卡断现象, 并且认为喉道长度及喉道宽度对卡断有一定影响.
相比于实验, 数值模拟方法可以更好的帮助理解及预测孔隙尺度气液两相流动状态及卡断现象. 在孔隙尺度上模拟多相流的方法有孔隙网络模型(PNM)[16]、光滑粒子动力学方法(SPH)[17]、密度泛函流体动力学方法(DFH)[18]、流体体积方法(VOF)[19]以及格子玻尔兹曼方法[20]等. 部分学者采用VOF方法研究了缩颈孔喉系统中的卡断现象, Raeini等[21]利用VOF方法研究了星形喉道截面的卡断现象, 并考虑了动态因素的影响. Starnoni和Pokrajac[22]采用VOF方法考虑接触角及黏度比对卡断影响, 研究表明, 喉道截面圆度越高发生卡断的临界接触角越小, 黏度比的增大也会降低发生卡断的临界接触角. Zhang等[23]采用VOF方法研究了方形截面喉道的卡断现象, 并认为高毛管数及高黏度比可以有效抑制卡断的发生, Cha等[13]利用VOF方法模拟了矩形截面的卡断现象, 并给出了卡断发生的孔喉宽度比条件. 然而, VOF方法很难准确计算界面的法向和曲率, 且界面重构方法复杂, 向高维推广困难[24]. 格子玻尔兹曼方法是介于微观分子动力学模拟方法和宏观传统数值模拟方法之间的一种介观数值模拟方法, 与VOF方法等界面跟踪或界面捕捉方法相比, 格子玻尔兹曼方法中的气液界面可以自动产生, 演化与迁移, 无须采用任何界面跟踪或界面捕捉技术[25], 并能直接处理流体与流体间的作用力及流体与固壁间的作用力[26]. 目前多相格子玻尔兹曼模型可以分为四类, 包括颜色梯度模型[27]、伪势模型[28]、自由能模型[29]以及相场模型[30]. 张磊等[31]基于颜色梯度模型在孔喉系统中模拟非混相驱替过程, 并认为多孔介质的非均质性越强, 流体越容易发生卡断. 赵玉龙等[32]基于颜色梯度模型, 模拟地层高温高压条件下致密气驱替地层水的流动过程, 并发现在多孔介质中大量连通的微小通道被卡断水占据. Alpak等[33]采用相场模型模拟3D缩颈孔隙中的两相驱替过程, 观察到了卡断现象并与Roof[6]静态准则保持较好的一致性, 但是模拟两相流体的黏度相差不大, 不能反应气液两相的黏度比. Wei等[34]采用伪势模型模拟单个气泡通过喉道过程, 并观察到了卡断现象. Zhang等[35]采用改进的双组分伪势模型模拟气驱水过程, 成功在喉道壁面捕捉到了水膜的存在, 并认为微纳尺度气水两相流动过程中壁面液膜的厚度不可忽略.
但是, 目前对卡断现象的研究往往建立在缩颈渐变孔喉系统中, 弱化了壁面液膜厚度的影响, 同时宏观VOF方法不能充分表征壁面和液膜相互作用. 为此, 本文在原始的伪势格子玻尔兹曼方法的基础上, 改进了流体之间作用力格式并添加流固作用力以捕捉气液两相驱替过程中喉道壁面上的液膜变化规律及其对气液两相流动状态的影响. 并进一步研究了驱替压差(毛管数Ca), 喉道长度及宽度等因素对液膜厚度及流动状态的影响规律. 本文的研究内容为表征复杂多孔介质的卡断机理提供了理论及模拟基础.
1. 单组分多相伪势模型
1.1 经典BGK碰撞模型
在格子玻尔兹曼方法中, 采用离散的密度分布函数表征流体的运动规律, 标准BGK碰撞过程可以表示为[36]
$$ \begin{split} & {f_i}({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t + \Delta t) - {f_i}({\boldsymbol{x}},t) = \hfill \\ &\qquad {\text{ }} - \frac{{\Delta t}}{\tau }[{f_i}({\boldsymbol{x}},t) - f_i^{\text{eq}}({\boldsymbol{x}},t)] + \Delta {f_i}({\boldsymbol{x}},t) \end{split} $$ (1) 式中, fi(x,t)为i方向密度分布函数, Δx和Δt分别表示网格步长和时间步长, τ表示松弛时间, 其取值与流体的运动黏度ν相关, 表示为
$$ \tau = \frac{\nu }{{c_s^2\Delta t}} + 0.5 $$ (2) ei为各方向的离散速度, 对于D2Q9模型表示为[37]
$$ {{\boldsymbol{e}}_i} = \left\{ \begin{gathered} 0,{\text{ }}i = 0 \hfill \\ \left\{ {\cos \left[ {\frac{{(i - 1){\text{π}} }}{2}} \right],\sin \left[ {\frac{{(i - 1){\text{π}} }}{2}} \right]} \right\},{\text{ }}i = 1\sim 4 \hfill \\ \sqrt 2 \left\{ {\cos \left[ {\frac{{(i - 5){\text{π}} }}{2} + \frac{{\text{π}} }{4}} \right],\sin \left[ {\frac{{(i - 5){\text{π}} }}{2} + \frac{{\text{π}} }{4}} \right]} \right\},{\text{ }}i = 5\sim 8 \hfill \\ \end{gathered} \right. $$ (3) 平衡密度分布函数
$f_i^{\text{eq}}({\boldsymbol{x}},t)$ 采用下式计算$$ f_i^{{\text{eq}}}({\boldsymbol{x}},t) = {w_i}\rho \left[1 + \frac{{{{\boldsymbol{e}}_i} \cdot {\boldsymbol{u}}}}{{c_s^2}} + \frac{{{{({{\boldsymbol{e}}_i} \cdot {\boldsymbol{u}})}^2}}}{{2c_s^4}} - \frac{{{{\boldsymbol{u}}^2}}}{{2c_s^2}}\right] $$ (4) 式中, wi为权重因子, 对于D2Q9模型, i = 0, wi = 4/9, i = 1 ~ 4, wi = 1/9, i = 5 ~ 8, wi = 1/36; cs为当地声速, 数值上
${c_s} = 1/\sqrt 3 $ 宏观密度ρ和速度u的计算如下$$ \rho = \sum\limits_i {{f_i}}\;\;\;\; $$ (5) $$ {\boldsymbol{u}} = \frac{{\displaystyle\sum\limits_i {{{\boldsymbol{e}}_i}{f_i}} }}{\rho } $$ (6) 1.2 流体之间作用力及流固作用力
原始伪势模型中流体之间的作用力定义为[28]
$$ {{\boldsymbol{F}}_{\rm{int} }}({\boldsymbol{x}},{\boldsymbol{x'}}) = - G(\left| {{\boldsymbol{x}} - {\boldsymbol{x'}}} \right|){\psi _\sigma }({\boldsymbol{x}}){\psi _{\bar \sigma }}({\boldsymbol{x}})({\boldsymbol{x'}} - {\boldsymbol{x}}) $$ (7) 式中, G为格林函数, ψ(x,t)为有效质量, 其定义为ψ(ρ) = ρ0[1−exp(−ρ/ρ0)], ρ0 = 1为参考密度. 进一步, 考虑第一层粒子之间的作用力, 上述作用力格式可以简化为
$$ {{\boldsymbol{F}}_{\rm{int} }}({\boldsymbol{x}},t) = - g\psi ({\boldsymbol{x}},t)c_s^2\sum\limits_i {w({{\left| {{{\boldsymbol{e}}_i}} \right|}^2})} \psi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t){{\boldsymbol{e}}_i} $$ (8) 式中, g为流体之间作用强度的控制参数, w(|ei|2)为权重系数, |ei|2 = 1, w(|ei|2) = 1/3; |ei|2 = 2, w(|ei|2) = 1/12. 采用原始的伪势模型的作用力格式会在两相界面处存在较大的虚假速度. 因此, 本文采用Gong等[38]提出的作用力改进格式以降低虚假速度, 具体表示为
$$ {{\boldsymbol{F}}_{\rm{int} }} = - \beta g\psi \nabla \psi ({\boldsymbol{x}}) - \frac{{1 - \beta }}{2}g{\nabla ^2}\psi ({\boldsymbol{x}}) $$ (9) 进一步考虑最近及次近的节点之间作用力, Chen等[39]提出了上述作用力的简化格式
$$ \begin{split} {{\boldsymbol{F}}_{\rm{int} }} = & - \frac{{1 - \beta }}{2}gc_s^2\sum\limits_{i = 1}^N {w({{\left| {{{\boldsymbol{e}}_i}} \right|}^2})} {\psi ^2}({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t){{\boldsymbol{e}}_i} - \hfill \\ & {\text{ }} \beta \frac{g}{2}\psi ({\boldsymbol{x}})c_s^2\sum\limits_{i = 1}^N {w({{\left| {{{\boldsymbol{e}}_i}} \right|}^2})\psi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t){{\boldsymbol{e}}_i}} \hfill \\ \end{split} $$ (10) 式中, β为调节参数, 通过改变β大小可以改善模型的热力学一致性并降低虚假速度. 若仅考虑第一层节点的作用力, 可将式(9)进一步扩展为[40]
$$\begin{split} {{\boldsymbol{F}}_x} =& gc_s^2\left\{ { - \beta {\psi _{i,j}}\Big[{\gamma _1}({\psi _{i + 1,j}} - {\psi _{i - 1,j}}) + {\gamma _2}\Big({\psi _{i + 1,j + 1}} - {\psi _{i - 1,j + 1}}} \right. + \\ & {\psi _{i + 1,j - 1}} - {\psi _{i - 1,j - 1}}\Big)\Big] - \frac{{1 - \beta }}{2}\Big[{\gamma _1}\Big(\psi _{i + 1,j}^2 - \psi _{i - 1,j}^2\Big) + \\ & {\gamma _2}\Big(\psi _{i + 1,j + 1}^2 - \psi _{i - 1,j + 1}^2 + \psi _{i + 1,j - 1}^2 - \psi _{i - 1,j - 1}^2\Big)\Big]\Big\}\\[-12pt] \end{split}$$ (11) $$\begin{split} {{\boldsymbol{F}}_y} =& gc_s^2 \Big\{ - \beta {\psi _{i,j}}\Big[{\gamma _1}({\psi _{i,j + 1}} - {\psi _{i,j - 1}}) + {\gamma _2}\Big({\psi _{i + 1,j + 1}} - {\psi _{i + 1,j - 1}} +\hfill \\ & {\psi _{i - 1,j + 1}} - {\psi _{i - 1,j - 1}}\Big)\Big] - \frac{{1 - \beta }}{2}\Big[{\gamma _1}(\psi _{i,j + 1}^2 - \psi _{i,j - 1}^2) +\hfill \\ & {\gamma _2}\Big(\psi _{i + 1,j + 1}^2 - \psi _{i + 1,j - 1}^2 + \psi _{i - 1,j + 1}^2 - \psi _{i - 1,j - 1}^2\Big)\Big] \Big\}\\[-12pt]\end{split} $$ (12) 式中, γ1和γ2为作用力格式的离散系数, 其中γ1 = 1/3, γ2 = 1/12. 上述离散作用力格式也称为E4作用力格式. 读者可以参考文献[40]的研究, 将作用力扩展至E6或E8作用力格式获得更低的虚假速度及更精确的界面张力.
流体与壁面的作用力类似于流体之间的作用力格式, 其定义如下[41]
$$ {{\boldsymbol{F}}_{{\text{ads}}}}({\boldsymbol{x}},t) = - g\psi ({\boldsymbol{x}},t)\sum\limits_i {{w_i}} \psi ({\rho _w})s({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t){{\boldsymbol{e}}_i} $$ (13) 式中, s(x + eiΔt,t)为判断流体节点和壁面的标识函数, 如果x + eiΔt是固体节点其值为1, x + eiΔt是流体节点则其值为0. ρw为壁面密度, 通过调节ρw的大小可以模拟不同的接触角.
1.3 状态方程
原始的伪势模型的状态方程可以通过泰勒展开获得[28], 具体表示为
$$ p = c_s^2\rho + \frac{{c_s^2g}}{2}{\psi ^2} $$ (14) 当采用上述的状态方程时, 会导致严重的热力学不一致性, 同时模拟的最大密度比和黏度比均在10以内. 文献[42]提出将真实的流体状态方程与伪势模型结合, 大大提高了伪势模型的热力学一致性. 本文采用常用的Redlich-Kwong (RK)状态方程表征气液两相, RK状态方程的定义为
$$ {p_{{\text{EOS}}}} = \frac{{\rho RT}}{{1 - b\rho }} - \frac{{a{\rho ^2}}}{{\sqrt T (1 + b\rho )}} $$ (15) 式中, pEOS为实际流体的压力, R为气体常数, T为系统的温度, ρ为实际流体的密度, a, b为临界参数, 其中a = 0.42748R2Tc2.5/pc, b = 0.08662RTc/pc. Huang等[43]给出了伪势模型中的参数取值, a = 2/49, b = 2/21, R = 1, 临界温度Tc = 0.196 1, 临界压力pc = 0.178 4, 临界密度ρc = 2.988 7. 在添加RK状态方程后可以采用下述方法计算伪势
$$ \psi = \sqrt {\dfrac{{2\Biggr[\dfrac{{\rho RT}}{{1 - b\rho }} - \dfrac{{a{\rho ^2}}}{{\sqrt T \Big(1 + b\rho \Big)}} - c_s^2\rho \Biggr]}}{{gc_s^2}}} $$ (16) 值得注意的是, 当采用上述格式计算伪势时, g将会在计算压力p的时候消掉, 不再具有实际物理意义, 仅作为保证根号内的符号为正的作用, 故本文的模拟中取g = −1.
1.4 外力项处理格式
在添加了流体-流体作用力及流固作用力后, 本文采用Kupershtokh等[44]提出的精确差分方法(EDM)格式处理外力项, 其定义如下
$$ \Delta {f_i}({\boldsymbol{x}},t) = f_i^{\text{eq}}\Biggr(\rho ,{\boldsymbol{u}} + \frac{{{{\boldsymbol{F}}_{{\text{total}}}}\Delta t}}{\rho }\Biggr) - f_i^{{\text{eq}}}(\rho ,{\boldsymbol{u}}) $$ (17) 式中, Ftotal表示流体之间作用力与流固作用力之和, EDM格式能够适应较宽的温度范围, 是目前最常见的外力处理格式之一.
1.5 边界条件
格子玻尔兹曼方法中的边界条件对数值计算精度, 计算稳定性以及计算效率具有很大影响[45]. 常用于模拟两相流动边界条件包括周期边界以及无滑移反弹边界等. 其中, 周期边界指流体粒子从一侧边界离开流场时, 在下一时间步会从流场的另一侧流入流场, 周期边界可以保证整个系统的质量和动量守恒. 传统的周期边界采用增加虚拟流体节点(x0, xN + 1)的方法实现, 可以表示为[46]
$$ \left. \begin{array}{l} {f_{1,5,8}}({x_0},y,t) = {f_{1,5,8}}({x_N},y,t) \hfill \\ {f_{3,6,7}}({x_{N + 1}},y,t) = {f_{3,6,7}}({x_1},y,t) \end{array} \right\} $$ (18) 对于静止的固体边界, 常用的处理方法就是对边界上的粒子进行反弹处理, 反弹边界可以严格保证边界上速度无滑移, 反弹边界可以表示为[46]
$$ \left. \begin{array}{l} {f_{2,5,6}}(x,{y_{{\text{bottom}}}},t + \Delta t) = {f_{4,7,8}}(x,{y_{{\text{bottom}}}},t) \hfill \\ {f_{4,7,8}}(x,{y_{{\text{top}}}},t + \Delta t) = {f_{2,5,6}}(x,{y_{{\text{top}}}},t) \end{array} \right\} $$ (19) 式中, ybottom和ytop分别为位于底部及顶部的固体边界.
2. 模型验证
2.1 热力学一致性验证
热力学不一致性一直是伪势模型中的主要问题之一, 将直接影响气液密度比及界面张力的计算精度. 因此, 需要首先评价模型的热力学一致性. 考虑到气液弯液面引起的毛管压力会对气液密度有一定影响, 采用假想的平直界面气液模型来研究气液共存密度, 并与麦克斯韦理论计算结果进行对比来评价模型的热力学一致性.
首先, 建立31 × 201的格子空间, 在x和y方向均设置为周期边界, 只改变温度, 以模拟不同温度下的气液共存密度, 采用Huang等[47]密度初始化方法提高数值稳定性, 该方法可以表示为
$$ \begin{split} \rho (y) =& {\rho _v} + \frac{{{\rho _l} - {\rho _v}}}{2}{\text{abs}}\left[ \tanh \frac{{2(y - 50)}}{W} -\right. \hfill \\ & \left. \tanh \frac{{2(y - 150)}}{W} \right] \end{split} $$ (20) 式中, W为相界面宽度, 一般为2 ~ 5个格子, ρv和ρl分别为麦克斯韦理论的气相和液相密度. 采用上述密度初始化方法之后, 在50 ≤ y ≤150的位置为液相, 其他位置为气相, 气液界面为直线, 可以忽略毛管压力作用[47]. 对比β = 1和β = 1.125及β = 1.25三种不同参数条件下的气液共存密度和麦克斯韦理论计算的气液共存密度, 结果如图1所示, β = 1.125时, LB模型与麦克斯韦理论所得的气液共存密度基本一致, 据此采用β = 1.125.
2.2 界面张力验证
在验证热力学一致性之后, 通过模拟静态圆形液滴来计算界面张力. 首先, 构建3种不同网格尺寸(101 × 101, 151 × 151, 201 × 201)的格子空间, 模拟静态圆形液滴密度分布并进行网格独立性检验. 采用Huang等[48]提出的密度初始化方法
$$ \begin{split} \rho (x,y) = &\frac{{{\rho _l} + {\rho _v}}}{2} -\hfill \\ & \frac{{{\rho _l} - {\rho _v}}}{2} \tanh \left\{\frac{{2\left[\sqrt {{{(x - {x_0})}^2} + {{(y - {y_0})}^2}} - {R_0}\right]}}{W}\right\} \end{split} $$ (21) 式中, R0为液滴初始化半径, (x0, y0)为液滴初始圆心位置, 模拟条件选取为: T = 0.8Tc, ρv = 0.342, ρl = 6.60, 四周采用周期边界, 其余参数与2.1节中保持一致. 为了便于比较, 将不同网格尺寸中液滴的直径与网格尺寸的比值保持一致, 模拟时间步选为50000步, 足以使模拟结果达到稳定. 模拟结果如图2所示, 不同网格尺寸下中心线处(x = Nx/2)对应的气相及液相密度分布几乎一致, 说明模拟结果具有网格独立性. 模拟的密度分布结果仅在两相过渡区域存在微小偏差, 当网格密度超过151 × 151时, 密度分布的模拟结果几乎一致, 足以满足模拟精度要求, 为了节约计算资源量, 选取网格尺寸为151 × 151的格子空间进行界面张力模拟验证.
进一步模拟T = 0.7Tc, T = 0.75Tc以及T = 0.8Tc三种不同温度条件下圆形液滴的静态行为, 模拟达到平衡后根据等式(14)计算实际的气液压差, 利用Young-Laplace方程∆P = Pin−Pout = σ/R确定界面张力[49]. 本文认为气液的分界面处的密度为(ρv + ρl)/2, 并以此为根据确定液滴直径大小, 通过改变液滴的不同直径, 获得一系列的压差和液滴半径的对应关系. 如图3所示, 压差与半径的倒数呈现良好的线性关系(R2均大于0.999), 通过线性拟合得到不同温度下气液界面张力大小分别为0.463 lu, 0.361 lu, 0.264 lu (其中lu代表格子单位).
2.3 静态接触角验证
本节在考虑流体-流体的作用力的基础上, 添加流体与固体壁面之间的作用力, 模拟液滴与壁面间的不同接触角. Chen等[39]给出了一种液滴在壁面上的接触角计算方法
$$ r = \frac{{4\xi _2^2 + \xi _1^2}}{{8{\xi _2}}} $$ (22) $$ \theta = \left\{ \begin{gathered} \arcsin \left(\frac{{{\xi _1}}}{{2r}}\right),{\text{ }}\theta \leqslant {90^ \circ } \hfill \\ {\text{π}} - \arcsin \left(\frac{{{\xi _1}}}{{2r}}\right),{\text{ }}\theta > {90^ \circ } \hfill \\ \end{gathered} \right. $$ (23) 式中, θ为液滴的润湿角, ξ1为液滴与固体壁面接触线的长度, ξ2为液滴的高度. 为了模拟验证不同的润湿角, 建立151 × 151的格子区域, 在x方向设置周期边界, 在y = 0及y = Ny设置为固体壁面, 并采用反弹边界[50], 模拟温度T=0.8Tc, 密度初始化方法采用Li等[51]提出的密度初始化方法, 其具体格式如下
$$ \rho (x,y)=\left\{\begin{array}{l}{\rho }_{l},\;{\rm{if }}\;\;{(x - {x_0})^2} + {(y - {y_0})^2}\leqslant {R}_{0}^{2}\\ {\rho }_{v}\text{, }\text{otherwise}\end{array} \right.$$ (24) 式中, R0为液滴初始化半径, 本文取15 lu, x0, y0为液滴初始化圆心位置, 本文设置(x0, y0) = (75, 135). 模拟结果如图4所示, 通过改变不同的壁面密度ρw可以模拟0° ~ 180°的润湿角.
2.4 角隅液相滞留特征
在方形截面孔隙及喉道中, 气相驱替液相之后, 液相将会在孔隙及喉道的角隅处滞留, 液相的滞留特征(曲率半径)对非圆形孔隙中的卡断机理有一定的影响[52]. 为此, 本文采用提出的伪势模型模拟液相在角隅的曲率半径与液相饱和度Sl的关系, 验证模型的准确性. 建立151 × 151的格子空间, 模拟温度T = 0.8Tc, 四周均为固体壁面并采用反弹边界, 润湿角设置为0°, 其余参数与2.3节中一致. 在文献[51]基础上提出一种新的初始化密度方法以表征液相在角隅中的不同曲率半径, 具体格式为
$$ \rho (x,y) = \left\{ \begin{gathered} {\rho _l},{\text{ if }}{(x - {x_0})^2} + {(y - {y_0})^2} \geqslant R_1^2 \hfill \\ {\rho _v},{\text{ otherwise}} \hfill \\ \end{gathered} \right. $$ (25) 式中, (x0,y0) = (75,75)为气泡初始圆心坐标位置, R1为初始气泡半径, 将R1的取值范围设为75 ~
${\rm{75}}\sqrt 2 $ 之间即可表征角隅液的不同曲率半径(0 ~ 75). 如图5所示, 所提模型可以捕捉到壁面的液膜及液相在角隅处的曲率半径. 壁面上静态液膜厚度的表征方法可以参考文献[53-57], 本文主要聚焦评价液相在角隅处的不同曲率直径d及饱和度的Sl关系.通过将液相在角隅的曲率半径与液相饱和度模拟结果与Li等[58]提出的方形孔隙角隅液的曲率半径与液相饱和度的关系相对比来验证模型的准确性, 若不考虑壁面液膜厚度的影响, 理论模型可以简化为
$$ {S_l} = \left(1 - \frac{{\text{π}} }{4}\right){\varepsilon ^2} $$ (26) 式中, ε为无量纲曲率直径, ε = d/(L1−2h), L1为方形孔隙的长度, h为液膜厚度, d为角隅液的曲率直径. 图6给出了模拟结果与理论结果的对比, 可以看出, 模拟结果与理论结果基本一致, 说明所提模型可以准确评价液相在角隅的赋存规律以及确定角隅液的曲率半径.
![]()
图 6 角隅液曲率直径与饱和度的关系Figure 6. Relationship between corner liquid curvature diameter and saturation2.5 网格独立性检验
通过模拟不同网格尺寸条件下单管内气相驱替液相过程(弯液面顶点位置随时间变化), 对模型进行网格独立性检验. 建立4种不同的网格尺寸(241 × 41, 301 × 51, 361 × 61, 421 × 71), 将y方向的底部及顶部设置为固体边界并采用无滑移反弹格式. 在入口和出口分别设置为5~10格子的气相缓冲区域, 其余位置设置为液相, 在每个格子节点添加外力驱动流动, 在流动方向(x方向)采用周期边界, 当液相被驱替至出口的气相缓冲区域中时, 液相密度在下一个时间步被自动替换为气相密度, 实现缓冲区域内液相“消失”, 从而实现驱替过程[35]. 模拟条件选取为: T = 0.8Tc, ρv = 0.342, ρl = 6.60, 将模拟的驱替压差∆p以及模拟时间步长与取点时间间隔的比值保持一致. 模拟结果如图7所示, 对于不同的网格尺寸, 模拟结果几乎一致. 因此, 在后续驱替的模拟中, 选取的模拟条件网格尺寸为301 × 51.
3. 动态卡断机理及影响因素分析
3.1 静态卡断准则评价
Ransohoff等[8]基于Roof提出的卡断判断机制, (即喉道毛管力Pc-neck大于前缘毛管力Pc-front), 提出了方形孔喉截面的静态卡断准则. 其核心假设是在方形截面孔喉系统中, 液相(润湿相)在角隅处以角流为主要流动方式, 如图8所示, 忽略壁面液膜的影响, 缩颈方形截面孔喉系统的卡断判断准则可以表示为
$$ {P_{\text{c-neck}}} - {P_{\text{c-front}}} = \sigma \left(\frac{1}{{{R_t}}} + \frac{1}{{{R_{zt}}}} - \frac{{{C_m}}}{{{R_c}}}\right) > 0 $$ (27) 式中, Rc为孔隙曲率半径, Rt为喉道曲率半径, Rzt为喉道横向曲率半径, Cm为无量纲曲率, 与孔隙截面形状有关. 基于最小表面能模型, Ransohoff等[8]给出了方形孔隙中无量纲曲率Cm=1.89. 当孔隙收缩性较平缓时1/Rzt可以忽略, 式(27)进一步简化为Rc > 1.89Rt.
Ransohoff等[8]提出的静态准则一般适用于驱替压差较小的准静态情况. 为了验证该静态准则的适用性, 本文采用改进的伪势模型进行模拟验证. 为了方便将模拟结果与Ransohoff静态准则进行对比, 首先定义无量纲孔喉长度比L*及孔喉宽度比R*表征孔隙结构的变化, L*和R*可以表示为
$$ {L^*} = L/{L_p}\tag{28a} $$ $$ {R^*} = {R_t}/{R_c} \tag{28b}$$ 式中, Lp为整个孔喉系统的长度, L为喉道长度. 基于此, 可以得到卡断静态准则为R*≤0.53.
建立301 × 51的格子空间, 在120 < x < 180及10 < y < 40的位置处设置宽度为30 lu的喉道(L* = 0.2, R* = 0.6). x方向的边界设置与2.5节中一致, 在y = 0及y = 50及喉道壁面设置为固体壁面, 润湿角设置为0°, 并采用反弹边界. 模拟的两相流体性质为: 气相密度ρv = 0.53 lu、液相密度ρl = 6.08 lu、气相黏度μv = 0.088 lu、液相黏度μl = 1.013 lu、两相界面张力为σ = 0.178 lu. 选取的模拟密度比(ρr = 11.5)与15 MPa, 350 K温度条件下的水(980 kg/m3)与甲烷(90 kg/m3)的密度比(ρr = 11)接近(NIST化学数据库).
为了更清楚的表示驱替过程, 笔者将整个驱替过程分为三个阶段(以时间为分界线). 第一阶段, 气相在左端的孔隙中的驱替阶段; 第二阶段, 气相侵入喉道的驱替阶段; 第三阶段, 气相突破喉道进入右端孔隙的驱替阶段. 在这三个阶段中, 选取气相饱和度Sg以及滞留液相饱和度Srl作为量化表征液相变化的参数. 气相饱和度Sg是指在驱替过程中气相体积分数与整个区域体积分数之比; 滞留液相饱和度Srl是指在驱替过程中, 以两相弯液面的顶点为分界线(垂直线), 将整个区域分为两部分, 分界线以内的液相体积分数与整个区域体积分数的比值. 换句话说, Sg + Srl等效为石油工程领域中活塞式驱替(界面为平直界面)过程中气驱波及的全部面积. 捕捉了驱替过程中的Sg以及Sg + Srl随时间的变化规律以及三种阶段某一时刻下的两相流体的密度分布(图9标注). 第一阶段(t < 9000), 气相在左端的孔隙中的驱替阶段, 驱替过程中形成相对稳定的弯液面, Srl保持较低的数值; 第二阶段(9000≤ t≤13000), 气相侵入喉道的驱替阶段, 此时在孔隙角隅处及喉道壁面上存在滞留的液相, 随着模拟时间步的增大Srl逐渐增大; 第三阶段(t > 13000), 气相突破喉道进入右端孔隙的驱替阶段, 此时由于壁面的强润湿性, 气相无法与壁面接触, 导致右端孔隙壁面上滞留大量液相, Srl迅速增大, 并且孔隙右侧滞留的液相(润湿相)随时间逐渐发生回流到喉道中, 最终卡断相界面. 此外, 所提模型孔喉条件为(R* = 0.6), 并不满足卡断发生的静态准则条件(R*≤ 0.53), 但是通过模拟可以观察到卡断现象. 这是因为在驱替的第二和第三阶段, 气相无法与喉道壁面以及右端孔隙壁面接触, 导致气相突破喉道后的前缘曲率半径始终小于孔隙半径, 使得 Cm取值偏大, 从而导致基于角流推导的静态准则预测结果偏保守.
进一步计算了驱替过程中的第三阶段内的Pc-neck 及Pc-front数值变化. 如图10所示, 随着时间步的增大Pc-neck不断增大, Pc-front不断减小, 这是因为Srl增大使得Pc-neck不断增大, 气相突破喉道之后前缘曲率半径不断增大导致Pc-front不断减小, 最终Pc-neck > Pc-front, 使得滞留的液相发生回流从而卡断相界面.
![]()
图 10 不同时间步下Pc-neck 及Pc-front随时间变化Figure 10. Pc-neck and Pc-front changes at different time steps上述结果与之前静态准则判断发生卡断的条件一致, 说明Roof[6]基于毛管压力不平衡的卡断判断机制在突变孔喉系统中同样适用. 同时, 上述2D模拟结果与Wu等[15]在3D孔喉结构中气泡卡断实验结果基本一致, 说明模拟结果对实际的3D情况仍然有效. 然而, 大量研究表明, 动态因素对卡断现象有较大影响[11], 但是目前的理论模型都是准静态几何准则, 无法考虑动态因素的影响. 为此, 利用提出的伪势模型, 采用数值模拟方法评价毛管数Ca, L*及R*对气液两相驱替中相界面卡断的影响规律.
3.2 毛管数的影响
毛管数Ca是表征两相流动的重要参数, 其定义为
$$ Ca = \frac{{{\rho _n}{\nu _n}{u_{{\text{ave}}}}}}{\sigma } $$ (29) 式中, ρn, νn分别代表非润湿相的密度及运动黏度, 在本文的模拟中与气相的密度及运动黏度相同, uave为两相驱替的平均速度[9].
通过改变驱替压差大小进而改变毛管数, 研究毛管数对于卡断的影响. 建立301 × 51的格子空间, 在100 < x < 200及10 < y < 40的位置处设置长度为100 lu, 宽度为30 lu的喉道(L* = 0.33, R* = 0.6), 驱替压差∆p = 0.069, 其余条件均与3.1中相同. 每隔500时间步, 捕捉相界面前缘位置, 通过相界面位置与时间关系曲线获得两相驱替平均速度(uave = ∆s/∆t). 需要指出的是, ∆s和∆t表示发生卡断现象之前的位移和时间. 图11(a)给出了Ca = 4.8 × 10−3时两相驱替过程中的Sg以及Sg + Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布, 第一阶段的模拟结果与3.1节中一致; 当气相进入喉道之后, 孔隙角隅及喉道壁面始终存在滞留的液相, Srl逐渐增大. 由于驱替压差较小, 气相不能突破喉道, 最终气液相界面在喉道内保持平衡不再延伸, Srl保持不变(Srl = 0.11). 造成上述现象的原因是气液界面的毛管压力以及固液界面的作用力形成较大的流动阻力, 这种现象称为毛管“钉扎”效应[59]. 因此, 发生卡断的首要条件即为驱替压差可以克服毛管“钉扎”效应, 形成有效驱替.
增大驱替压差至∆p = 0.078, 图11(b)给出了Ca = 5.6 × 10−3时两相驱替过程的Sg以及Sg + Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布. 驱替的第一阶段与前文模拟结果一致; 在气相进入喉道之后, 孔隙角隅及喉道壁面始终存在滞留的液相, Srl逐渐增大, 且滞留在喉道壁面上的液相厚度呈非均匀分布. 随着模拟时间步的增大, 滞留的液相逐渐回流并卡断相界面. 值得一提的是, 发生卡断时间属于驱替的第二阶段, 说明在非渐变孔喉系统中卡断可以在气相突破喉道之前发生.
继续增大驱替压差至∆p = 0.105, 其余条件保持不变, 图11(c)给出了Ca = 8.0 × 10−3时两相驱替过程中的Sg以及Sg + Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布. 此时毛管数相对较大, 相比于之前的模拟结果, 在第二阶段, 喉道壁面上滞留的液相饱和度Srl较小; 随着模拟的时间步的增大, 气相不断向前流动, 在突破喉道前并未发生卡断, 直至气相到达出口端, 随着模拟时间步的继续增大, 滞留的液相发生回流并形成卡断. 并且随着毛管数的增大, 发生卡断的位置向喉道出口端移动, 这与Wei等[60]采用3D孔喉结构研究不同毛管数条件下气泡卡断的实验结果一致.
进一步增大驱替压差至∆p = 0.114, 其余条件保持不变, 图11(d)给出了Ca = 8.6 × 10−3时两相驱替过程中的Sg以及Sg + Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布. 此时驱替压差很大, 相比于之前的模拟结果Srl较小; 随着时间步的增大, 气相从喉道中突破, 并最终到达出口端. 在气相到达出口端以后, 随着模拟时间步的增大, 并未发生卡断现象. 这说明对于固定的孔喉结构, 存在一个临界毛管数, 当系统的毛管数大于临界毛管数时, 将会抑制滞留的液相在喉道内发生回流, 从而抑制卡断现象的发生.
图12给出了50000时间步内, 上述四种条件下相界面前缘顶点位置随时间的变化. 发生卡断之后, 把回退到喉道中的弯液面顶点作为实际的相界面位置, 这样当相界面位置突然变化时, 即可判断为发生卡断. 可以看出, 气相为连续相时, 驱替过程中可以发生连续卡断, 并且每次发生卡断的位置几乎不会变, 该过程类似于液滴/气泡的形成过程[61-62]. 同时, 对于特定的孔喉系统, 只有在满足一定的毛管数范围内才会发生卡断, 较大的毛管数会抑制驱替过程中滞留液相的回流从而抑制卡断的发生, 较小的毛管数无法克服毛管“钉扎”效应形成有效驱替. 模拟结果与文献[9]在渐变孔喉中研究气泡在压力梯度作用下的卡断行为所得出的结论一致.
3.3 孔喉长度比的影响
孔喉结构对卡断同样具有重要影响[9-10], 为此本节利用提出的伪势模型研究喉道长度对卡断的影响. 建立301 × 51的格子空间, 在10 < y < 40的位置处设置相同宽度(R* = 0.6), 不同长度的喉道来研究不同L*对卡断的影响. 对于不同的喉道长度, 驱替压差比毛管数更能体现动态因素对卡断的影响, 因此采用驱替压差代替毛管数表征动态参数的影响.
在模拟过程中, 改变驱替压差从0.036至0.123, 改变L*从0.08至0.4来确定卡断发生的临界条件; 这里仅用∆p = 0.084, L* = 0.3和∆p = 0.084, L* = 0.2两种模拟条件为例, 展示不同孔喉长度比的驱替效果. 图13给出了∆p = 0.084, L* = 0.3及L* = 0.2条件下, 不同时间步的驱替过程. 随着喉道长度的减小, 流动阻力减小, 两相驱替平均速度增大, 驱替过程中滞留在孔隙角隅及喉道壁面处的液相饱和度Srl较低. 值得一提的是, 设置的孔喉宽度比R* = 0.6, 不满足静态准则预测发生卡断的条件(R*≤0.53), 但是模拟结果可以观察到卡断现象, 说明即使不满足静态准则预测发生卡断的孔喉宽度比, 足够长的喉道长度也会促进卡断的发生. 因此, 在相同驱替压差的作用下, 不同喉道长度的孔喉系统会呈现不同的流动状态.
为了确定不同孔喉长度比的卡断发生条件, 开展了一系列模拟, 得到不同孔喉长度比下发生卡断的临界驱替压差, 如图14所示. 对于固定的喉道宽度的孔喉系统, 存在临界孔喉长度比(L* = 0.08), 使得无论驱替压差如何变化, 气液两相驱替过程中不会出现卡断现象, 这与Deng等[10]提出的短波长的渐变孔喉中卡断会被抑制的结论一致. 随着喉道长度增大, 不同的驱替压差会使得气液驱替过程中出现卡断, 连续流动以及无效驱替三种流动状态. 当驱替压差小于临界驱替压差下界时, 无法克服毛管“钉扎”效应, 形成无效驱替; 当驱替压差大于临界驱替压差上界时, 气相将保持连续流动的状态, 抑制卡断的发生.
![]()
图 14 不同L*条件下气液流动状态与驱替压差的关系Figure 14. The relationship between gas-liquid flow state and displacement pressure difference with different L*3.4 孔喉宽度比的影响
除了喉道长度, 喉道的宽度对于卡断现象也有一定影响, 本节通过固定L*改变R*, 研究不同喉道宽度对卡断的影响. 建立301 × 51的格子空间, 在125 < x < 175设置长度为50 lu的喉道(L* = 0.167), 其余参数设置与3.3节相同. 在模拟过程中, 改变驱替压差从0.033至0.1, 改变R*大小从0.52至0.68来确定卡断临界条件; 这里仅用∆p = 0.084, R* = 0.52和∆p = 0.084, R* = 0.6两组条件为例, 展示不同孔喉宽度比对气液两相驱替过程中的流动状态的影响. 图15给出了∆p = 0.084, R* = 0.52及R* = 0.6两种条件下的驱替过程. 随着喉道宽度的减小, 流动阻力显著增大, 两相驱替平均速度减小, 驱替过程中滞留在孔隙角隅及喉道壁面处的液相饱和度Srl较高. 值得一提的是, R* = 0.52满足静态准则预测发生卡断的条件(R*≤ 0.53), 模拟结果同样观察到了卡断现象, 验证了模拟结果的准确性. 对于R* = 0.6的情况下, 驱替过程中始终不会发生卡断现象, 这也说明了喉道宽度增大将会抑制卡断的发生.
为了确定不同孔喉宽度比对卡断现象的影响, 开展了一系列模拟, 得到了不同孔喉宽度比下发生卡断的临界驱替压差. 如图16所示, 与3.3节类似, 对于固定喉道长度的孔喉系统, 气液两相在不同驱替压差的作用下形成气相连续流动, 卡断及无效驱替3类流动状态. 随着喉道宽度的增加, 发生卡断的驱替压差范围逐渐减小, 当孔喉宽度比增大至R* = 0.68时, 无论施加多大的驱替压差, 喉道内都不会发生卡断现象. 值得一提的是, 对于R* = 0.52的孔喉结构(满足静态准则预测发生卡断的条件), 当驱替压差大于0.1时, 驱替过程中将不会发生卡断现象, 这说明即使达到静态准则所预测的卡断条件, 较大的驱压差也可以抑制卡断的发生. 同时, 基于等式(27)预测临界孔喉宽度比为(R* = 0.53), 相比于通过模拟得出的临界孔喉宽度比(R* = 0.68)降低28.3%.
![]()
图 16 不同R*条件下气液流动状态与驱替压差的关系Figure 16. The relationship between gas-liquid flow state and displacement pressure difference with different R*4. 结 论
本文在原始单组分伪势格子玻尔兹曼模型基础上改进流体-流体作用力格式, 建立了改进的伪势模型, 并在此基础上模拟了孔喉系统中的气液两相在不同驱替条件下的流动状态, 得出以下结论.
(1) 明确了非渐变孔喉结构卡断机理, 即孔喉毛管压力不平衡使得驱替过程中滞留的液相(润湿相)随时间逐渐发生回流是造成相界面卡断的根本原因. 但是, 模拟过程中在孔隙角隅及喉道壁面观察到大量滞留的液相, 这导致基于角流的假设思想建立的静态卡断判断准则将会高估喉道右侧气泡的曲率半径从而低估卡断发生的条件.
(2) 揭示了动态因素驱替压差(毛管数)对气液两相驱替的影响规律. 在非渐变孔喉系统中, 只有 驱替压差(毛管数)大小在一定范围时, 喉道内才会发生卡断; 超过该毛管数上界, 即使满足静态准则条件, 卡断也会被抑制; 低于该下界, 无法完成驱替; 同时毛管数的增大使得发生卡断的位置向喉道出口端移动.
(3) 揭示了孔喉长度比对气液两相驱替的影响规律. 对于固定宽度的孔喉系统, 即使不满足静态卡断发生的孔喉宽度比(R*≤ 0.53), 足够长的喉道也可以促进卡断的发生, 随着喉道长度增大, 发生卡断的驱替压差范围越大. 并且存在一个临界喉道长度, 使得无论驱替压差如何变化, 喉道内都不会发生卡断, 对于本文的模型, 临界孔喉长度比L* = 0.08.
(4) 揭示了孔喉宽度比对气液两相驱替的影响规律. 对于固定长度的孔喉系统, 喉道宽度越大, 发生卡断的驱替压差范围越小; 并且存在一个临界喉道宽度, 使得无论驱替压差如何变化, 喉道内都不会发生卡断. 对于本文的模型, 临界孔喉宽度比R* = 0.68. 若采用静态准则预测该模型发生卡断的条件(R* = 0.53), 将会低估28.3%.
-
![]()
图 6 角隅液曲率直径与饱和度的关系
Figure 6. Relationship between corner liquid curvature diameter and saturation
![]()
图 10 不同时间步下Pc-neck 及Pc-front随时间变化
Figure 10. Pc-neck and Pc-front changes at different time steps
![]()
图 14 不同L*条件下气液流动状态与驱替压差的关系
Figure 14. The relationship between gas-liquid flow state and displacement pressure difference with different L*
-
[1] Yun W, Kovscek AR. Microvisual investigation of polymer retention on the homogeneous pore network of a micromodel. Journal of Petroleum Science Engineering, 2015, 128: 115-127 doi: 10.1016/j.petrol.2015.02.004
[2] 柳占立, 庄茁, 孟庆国等. 页岩气高效开采的力学问题与挑战. 力学学报, 2017, 49(3): 507-516 (Liu Zhanli, Zhuang Zhuo, Meng Qingguo, et al. Problems and challenges of mechanics in shale gas efficient exploitation. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 507-516 (in Chinese) [3] 袁士义, 王强, 李军诗等. 注气提高采收率技术进展及前景展望. 石油学报, 2020, 41(12): 1623-1632 (Yuan Shiyi, Wang Qiang, Li Junshi, et al. Technology progress and prospects of enhanced oil recovery by gas injection. Acta Petrolei Sincia, 2020, 41(12): 1623-1632 (in Chinese) doi: 10.7623/syxb202012014 [4] 高云丛, 赵密福, 王建波等. 特低渗油藏CO2非混相驱生产特征与气窜规律. 石油勘探与开发, 2014, 41(1): 79-85 (Gao Yuncong, Zhao Mifu, Wang Jianbo, et al. Performance and gas breakthrough during CO2 immiscible flooding in ultra-low permeability reservoirs. Petroleum Exploration and Development, 2014, 41(1): 79-85 (in Chinese) [5] Kong D, Gao Y, Sarma H, et al. Experimental investigation of immiscible water-alternating-gas injection in ultra-high water-cut stage reservoir. Advances in Geo-Energy Research, 2021, 5(2): 139-152 doi: 10.46690/ager.2021.02.04
[6] Roof JG. Snap-off of oil droplets in water-wet pores. Society of Petroleum Engineers Journal, 1970, 10(1): 85-90 doi: 10.2118/2504-PA
[7] Gauglitz PA, StLaurent CM, Radke CJ. Experimental determination of gas-bubble breakup in a constricted cylindrical capillary. Industrial & Engineering Chemistry Research, 1988, 27(7): 1282-1291
[8] Ransohoff TC, Gauglitz PA, Radke CJ. Snap-off of gas bubbles in smoothly constricted noncircular capillaries. AIChE Journal, 1987, 33(5): 753-765 doi: 10.1002/aic.690330508
[9] Tsai TM, Miksis MJ. Dynamics of a drop in a constricted capillary tube. Journal of Fluid Mechanics, 2016, 274(274): 197-217
[10] Deng W, Cardenas MB, Bennett PC. Extended Roof snap-off for a continuous nonwetting fluid and an example case for supercritical CO2. Advances in Water Resources, 2014, 64: 34-46 doi: 10.1016/j.advwatres.2013.12.001
[11] Deng W, Balhoff M, Cardenas MB. Influence of dynamic factors on nonwetting fluid snap-off in pores. Water Resources Research, 2015, 51(11): 9182-9189 doi: 10.1002/2015WR017261
[12] Tian J, Kang YL, Xi ZY, et al. Real-time visualization and investigation of dynamic gas snap-off mechanisms in 2-D micro channels. Fuel, 2020, 279: 118232 doi: 10.1016/j.fuel.2020.118232
[13] Cha LM, Xie CY, Feng QH, et al. Geometric criteria for the snap-off of a non-wetting droplet in pore-throat channels with rectangular cross-sections. Water Resources Research, 2021, 57(7): e2020WR029476
[14] Tetteh JT, Cudjoe SE, Aryana SA, et al. Investigation into fluid-fluid interaction phenomena during low salinity waterflooding using a reservoir-on-a-chip microfluidic model. Journal of Petroleum Science and Engineering, 2021, 196: 108074 doi: 10.1016/j.petrol.2020.108074
[15] Wu YN, Fang SS, Dai CL, et al. Investigation on bubble snap-off in 3-D pore-throat micro-structures. Journal of Industrial and Engineering Chemistry, 2017, 54: 69-74 doi: 10.1016/j.jiec.2017.05.019
[16] Xiong QR, Baychev TG, Jivkov AP. Review of pore network modelling of porous media: Experimental characterisations, network constructions and applications to reactive transport. Journal of Contaminant Hydrology, 2016, 192: 101-117 doi: 10.1016/j.jconhyd.2016.07.002
[17] Monaghan JJ. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 1992, 30(1): 543-574 doi: 10.1146/annurev.aa.30.090192.002551
[18] Armstrong RT, Berg S, Dinariev O, et al. Modeling of pore-scale two-phase phenomena using density functional hydrodynamics. Transport in Porous Media, 2016, 112(3): 577-607 doi: 10.1007/s11242-016-0660-8
[19] Raeini AQ, Blunt MJ, Bijeljic B. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. Journal of Computational Physics, 2012, 231(17): 5653-5668 doi: 10.1016/j.jcp.2012.04.011
[20] Chen SY, Doolen GD. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics, 1998, 30(1): 329-364 doi: 10.1146/annurev.fluid.30.1.329
[21] Raeini AQ, Bijeljic B, Blunt MJ. Numerical modelling of sub-pore scale events in two-phase flow through porous media. Transport in Porous Media, 2014, 101(2): 191-213 doi: 10.1007/s11242-013-0239-6
[22] Starnoni M, Pokrajac D. Numerical study of the effects of contact angle and viscosity ratio on the dynamics of snap-off through porous media. Advances in Water Resources, 2018, 111: 70-85 doi: 10.1016/j.advwatres.2017.10.030
[23] Zhang CW, Yuan ZY, Matsushita S, et al. Interpreting dynamics of snap-off in a constricted capillary from the energy dissipation principle. Physics of Fluids, 2021, 33(3): 032112 doi: 10.1063/5.0044756
[24] 张嫚嫚, 孙姣, 陈文义. 一种基于几何重构的Youngs-VOF耦合水平集追踪方法. 力学学报, 2019, 51(3): 775-786 (Zhang Manman, Sun Jiao, Chen Wenyi. An interface tracking method of coupled Youngs-VOF and level set based on geometric reconstruction. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(3): 775-786 (in Chinese) [25] 李庆, 余悦, 唐诗. 多相格子Boltzmann方法及其在相变传热中的应用. 科学通报, 2020, 65(17): 1677-1693 (Li Qing, Yu Yue, Tang Shi. Multiphase lattice Boltzmann method and its applications in phase-change heat transfer. Chinese Science Bulletin, 2020, 65(17): 1677-1693 (in Chinese) doi: 10.1360/TB-2019-0769 [26] 臧晨强, 娄钦. 复杂微通道内非混相驱替过程的格子 Boltzmann方法. 物理学报, 2017, 66(13): 154-162 (Zang Chenqiang, Lou Qin. Lattice Boltzmann simulation of immiscible displacement in the complex micro-channel. Acta Physica Sinica, 2017, 66(13): 154-162 (in Chinese) [27] Rothman DH, Keller JM. Immiscible cellular-automaton fluids. Journal of Statistical Physics, 1988, 52(3): 1119-1127
[28] Shan XW, Chen HD. Lattice Boltzmann model for simulating flows with multiple phases and components. Physical Review E, 1993, 47(3): 1815 doi: 10.1103/PhysRevE.47.1815
[29] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulation of nonideal fluids. Physical Review Letters, 1995, 75(5): 830 doi: 10.1103/PhysRevLett.75.830
[30] He XY, Shan XW, Doolen GD. Discrete Boltzmann equation model for nonideal gases. Physical Review E, 1998, 57(1): R13 doi: 10.1103/PhysRevE.57.R13
[31] 张磊, 康立新, 景文龙等. 基于孔隙-喉道双通道模型的油液两相流动形态分析. 中国石油大学学报(自然科学版), 2020, 44(5): 89-93 (Zhang Lei, Kang Lixin, Jing Wenlong, et al. Flow behavior analysis of oil-water two-phase flow in pore throat doublet model. Journal of China University of Petroleum (Edition of Natural Science) , 2020, 44(5): 89-93 (in Chinese) [32] 赵玉龙, 刘香禺, 张烈辉等. 致密砂岩气藏气液流动规律及储层干化作用机理. 天然气工业, 2020, 40(9): 70-79 (Zhao Yulong, Liu Xiangyu, Zhang Liehui, et al. Laws of gas and water flow and mechanism of reservoir drying in tight sandstone gas reservoirs. Natural Gas Industry B, 2020, 40(9): 70-79 (in Chinese) [33] Alpak FO, Zacharoudiou I, Berg S, et al. Direct simulation of pore-scale two-phase visco-capillary flow on large digital rock images using a phase-field lattice Boltzmann method on general-purpose graphics processing units. Computational Geosciences, 2019, 23(5): 849-880 doi: 10.1007/s10596-019-9818-0
[34] Wei B, Hou J, Sukop MC, et al. Flow behaviors of emulsions in constricted capillaries: a Lattice Boltzmann simulation study. Chemical Engineering Science, 2020, 227: 115925 doi: 10.1016/j.ces.2020.115925
[35] Zhang T, Javadpour F, Li J, et al. Pore-scale perspective of gas/water two-phase flow in shale. SPE Journal, 2021, 26(2): 828-846 doi: 10.2118/205019-PA
[36] Bhatnagar PL, Gross EP, Krook MA model for collision processes in gases. I Small amplitude processes in charged and neutral one-component systems. Physical Review, 1954, 94(3): 511
[37] Qian YH, d'Humières D, Lallemand P. Lattice BGK models for Navier-Stokes equation. EPL (Europhysics Letters)
, 1992, 17(6): 479 doi: 10.1209/0295-5075/17/6/001 [38] Gong S, Cheng P. Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows. Computers & Fluids, 2012, 53: 93-104
[39] Chen L, Kang QJ, Mu YT, et al. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer, 2014, 76: 210-236 doi: 10.1016/j.ijheatmasstransfer.2014.04.032
[40] Mukherjee A, Basu DN, Mondal PK. Algorithmic augmentation in the pseudopotential-based lattice Boltzmann method for simulating the pool boiling phenomenon with high-density ratio. Physical Review E, 2021, 103(5): 053302 doi: 10.1103/PhysRevE.103.053302
[41] Huang HB, Li ZT, Liu SS, et al. Shan-and-Chen-type multiphase lattice Boltzmann study of viscous coupling effects for two-phase flow in porous media. International Journal for Numerical Methods in Fluids, 2009, 61(3): 341-354 doi: 10.1002/fld.1972
[42] Yuan P, Schaefer L. Equations of state in a lattice Boltzmann model. Physics of Fluids, 2006, 18(4): 042101 doi: 10.1063/1.2187070
[43] Huang HB, Sukop M, Lu XY. Multiphase Lattice Boltzmann Methods: Theory and Application. John Wiley & Sons, 2015
[44] Kupershtokh AL, Medvedev DA, Karpov DI. On equations of state in a lattice Boltzmann method. Computers & Mathematics with Applications, 2009, 58(5): 965-974
[45] 何雅玲, 王勇, 李庆. 格子Boltzmann方法的理论及应用. 北京: 科学出版, 2009 He Yaling, Wang Yong, Li Qing. Theory and Application of Lattice Boltzmann Method. Beijing: Science Press, 2009 (in Chinese))
[46] Krüger T, Kusumaatmaja H, Kuzmin A, et al. The lattice Boltzmann method. Springer International Publishing, 2017, 10: 4-15
[47] Huang JW, Yin XL, Barrufet M, et al. Lattice Boltzmann simulation of phase equilibrium of methane in nanopores under effects of adsorption. Chemical Engineering Journal, 2021, 419: 129625 doi: 10.1016/j.cej.2021.129625
[48] Huang HB, Krafczyk M, Lu XY. Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 2011, 84(2): 046710
[49] 史冬岩, 王志凯, 张阿漫. 一种模拟气液两相流的格子波尔兹曼改进模型. 力学学报, 2014, 46(2): 224-233 (Shi Dongyan, Wang Zhikai, Zhang Aman. A novel lattice boltzmann model simulating gas-liquid two-phase flow. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(2): 224-233 (in Chinese) [50] 胡五龙, 刘国峰, 晏石林等. 土壤水分布的孔隙尺度格子玻尔兹曼模拟研究. 力学学报, 2021, 53(2): 568-579 (Hu Wulong, Liu Guofeng, Yan Shilin, et al. Pore-scale lattice Boltzmann modeling of soil water Distribution. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 568-579 (in Chinese) [51] Li Q, Luo KH, Kang QJ, et al. Contact angles in the pseudopotential lattice Boltzmann modeling of wetting. Physical Review E, 2014, 90(5): 053301 doi: 10.1103/PhysRevE.90.053301
[52] Kovscek AR, Radke CJ. Gas bubble snap-off under pressure-driven flow in constricted noncircular capillaries. Colloids and Surfaces A:Physicochemical and Engineering Aspects, 1996, 117(1-2): 55-76
[53] Li J, Li XF, Wang XZ, et al. Water distribution characteristic and effect on methane adsorption capacity in shale clay. International Journal of Coal Geology, 2016, 159: 135-154 doi: 10.1016/j.coal.2016.03.012
[54] Li J, Li XF, Wu KL, et al. Thickness and stability of water film confined inside nanoslits and nanocapillaries of shale and clay. International Journal of Coal Geology, 2017, 179: 253-268 doi: 10.1016/j.coal.2017.06.008
[55] Li J, Li XF, Wu KL, et al. Water sorption and distribution characteristics in clay and shale: effect of surface force. Energy & Fuels, 2016, 30(11): 8863-8874
[56] 李靖, 李相方, 王香增等. 页岩黏土孔隙含水饱和度分布及其对甲烷吸附的影响. 力学学报, 2016, 48(5): 1217-1228 (Li Jing, Li Xiangfang, Wang Xiangzeng, et al. Effect of water distribution on methane adsorption capacity in shale clay. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5): 1217-1228 (in Chinese) [57] 李靖, 李相方, 李莹莹等. 储层含水条件下致密砂岩/页岩无机质纳米孔隙气相渗透率模型. 力学学报, 2015, 47(6): 932-944 (Li Jing, Li Xiangfang, Li Yingying, et al. Model for gas transport in nanopores of shale and tight formation under reservoir condition. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(6): 932-944 (in Chinese) [58] Li J, Chen ZX, Wu KL, et al. Effect of water saturation on gas slippage in circular and angular pores. AIChE Journal, 2018, 64(9): 3529-3541 doi: 10.1002/aic.16196
[59] 宋付权, 胡箫, 朱根民等. 纳米阵列中气体驱替液体的流动特征. 力学学报, 2018, 50(3): 553-560 (Song Fuquan, Hu Xiao, Zhu Genmin, et al. The characteristics of water flow displaced by gas in nano arrays. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(3): 553-560 (in Chinese) [60] Wei B, Wang YY, Wen YB, et al. Bubble breakup dynamics and flow behaviors of a surface-functionalized nanocellulose based nanofluid stabilized foam in constricted microfluidic devices. Journal of Industrial and Engineering Chemistry, 2018, 68: 24-32 doi: 10.1016/j.jiec.2018.07.025
[61] Si T. Dynamic behavior of droplet formation in dripping mode of capillary flow focusing. Capillarity, 2021, 4(3): 45-49 doi: 10.46690/capi.2021.03.01
[62] Liu HH, Zhang YH. Droplet formation in microfluidic cross-junctions. Physics of Fluids, 2011, 23(8): 082101 doi: 10.1063/1.3615643
-
期刊类型引用(0)
其他类型引用(1)
下载:
计量
- 文章访问数: 1204
- HTML全文浏览量: 395
- PDF下载量: 182
- 被引次数: 1
