0 引言

混溶驱替指2种或2种以上可混溶流体的驱替，既可以是单相流体(海水与淡水、地表污水与地下水、含示踪剂的流体和纯流体)，也可以是两相流体(CO₂与原油等)，其广泛应用于水文学、化学、医学和石油工程等领域^[1-4]，而混相驱替主要应用于石油工业领域，且更侧重于2种不同相态的可混溶流体的驱替^[3-7]。由于多孔介质孔隙结构特征(孔隙非均质性、孔隙连通性、迂曲度等)和流体流动特征(被驱替流体与驱替流体的黏度比、分子扩散、弥散、驱替流量等)的差异，混溶驱替过程常常表现出黏性指进现象。因此，掌握影响混溶驱替的主控因素对于弄清溶质运移规律和扩大驱替流体的体积波及效率具有重要意义。

国内外学者^[8-16]利用连续模型^[8]、随机模型^[9]、孔隙网络模拟^[10-11]和Lattice-Boltzmann模拟^[12-14]等方式开展了大量混溶驱替研究：Mohiuddin等^[8]利用Navier-Stokes方程、连续性方程和对流扩散方程研究了重力对混溶驱替过程的影响，该方法很难考虑多孔介质的孔隙结构特征，仅能进行定性分析，不能作定量评价；Tartakovsky等^[15]用粒子表示流体，运用光滑粒子流体力学方法(SPH)求解Lan‐ gevin方程，实现了混溶过程中对流和扩散的有效分离。Ebrahimi等^[9]运用小波变换法粗化非均匀多孔介质来实现混溶驱替的升尺度模拟，发现渗透率各向异性和黏度比均是影响油田尺度混溶驱替效率的关键因素；Xia^[12]重点研究了孔隙结构、Le数和Re数等对混溶驱替过程的影响；Stevenson等^[10]和Sideiqui等^[16]利用孔隙网络模型开展了流体的随机性或确定性运动模拟，研究了非均质性、流度比等因素对混溶驱的影响，但均未考虑弥散的作用。

为加深对混溶驱替微观机理的认识，通过对网络模型压力场和浓度场的耦合求解，开展考虑弥散的混溶驱替动态网络模拟，定量分析孔隙非均质性、孔隙连通性、黏度比和驱替流量等对岩心尺度下混溶驱替的影响，以期弄清混溶驱替过程的主控因素。

1 控制方程

当圆形管束中发生单相流体流动时，体积流量服从Hagen-Poiseuille方程^[17]，即

$ {q_{ij}} = \frac{{{{10}^9}\pi {\kern 1pt} r_{ij}^4({p_i} - {p_j})}}{{8\mu {\kern 1pt} l}} $

(1)

式中：q_ij为管束流量，m³/s；p_i，p_j分别为由管束所连接的节点i，j的压力，MPa；r_ij为管束半径，m；μ为流体黏度，mPa·s；l为管束长度，m。

当圆形管束中发生不混溶驱替(油和水)时，不同流体间存在明显的分界面，若驱替流体占据的管束长度为X_ij，则被驱替流体占据的长度为l - X_ij，此时管束的体积流量为

$ {q_{ij}} = \frac{{{{10}^9}\pi {\kern 1pt} r_{ij}^4({p_i} - {p_j})}}{{8[{\mu _1}{X_{ij}} + (l - {X_{ij}}){\mu _2}]}} $

(2)

式中：μ₁，μ₂分别为驱替流体、被驱替流体的流体黏度，mPa·s；X_ij为驱替流体或混合带所占据的管束长度，m；。

当圆形管束中发生混溶驱替(海水与淡水、CO₂与油)时，管束中不存在完全由驱替流体所占据的部分，而变为驱替流体与被驱替流体的混合带，此时管束的体积流量为^[17]

$ {q_{ij}} = \frac{{{{10}^9}\pi {\kern 1pt} r_{ij}^4({p_i} - {p_j})}}{{8[{\mu _{\rm{e}}}{X_{ij}} + (l - {X_{ij}}){\mu _2}]}} $

(3)

式中：μ_e为混合带的流体黏度，mPa·s。

混合带的黏度μ_e根据流体混合的四次方根理论^[18]进行计算

$ {\mu _e} = {\left( {\frac{{{C_s}}}{{\mu _1^{0.25}}} + \frac{{1 - {C_s}}}{{\mu _2^{0.25}}}} \right)^{ - 4}} $

(4)

式中：C_s为驱替流体所占据的体积分数(体积百分浓度)，m³/m³。

根据节点流量守恒，流入节点的流量应等于流出节点的流量，即

$ \sum {{q_{ij}}} = 0 $

(5)

根据式(3)-(5)可列出所有节点对应的线性方程组，将线性方程组转化为矩阵的形式[式(6)]，再利用共轭梯度法对矩阵求解得到网络模型的压力分布

$ A{\kern 1pt} {\kern 1pt} P = 0 $

(6)

式中：A为传导率系数矩阵；P为模型各个节点压力组成的向量。

在确定网络模型的初始压力分布后，需要对网络模型的浓度场分布进行求解。驱替流体在单个孔隙的流动满足Taylor-Aris弥散^[19-20]且服从一维对流弥散方程^[21-22](ADE)

$ {\frac{{\partial {\kern 1pt} C}}{{\partial {\kern 1pt} t}} = D\frac{{{\partial ^2}C}}{{\partial {\kern 1pt} {x^2}}} - {v_{\rm{m}}}\frac{{\partial {\kern 1pt} C}}{{\partial {\kern 1pt} x}}} $

(7)

$ {D = {D_m} + \frac{{r_{ij}^2v_{\rm{m}}^2}}{{48{\kern 1pt} {D_{\rm{m}}}}}} $

(8)

式中：C为体积百分浓度，m³/m³；t为时间，s；D为Taylor-Aris弥散系数，m²/s；v_m为管束的平均流速，m/s；D_m为分子扩散系数，m²/s。

采用显式差分格式求解非线性方程，计算管束中的浓度分布。兼顾模型求解的准确性和计算的可行性，将单根管束剖分为4个网格(以二维网络模型为例，图 1)，则空间步长△x等于管束长度的1/4。

图 1中，若驱替流体从节点1(i，j)流入节点5 (i+1，j)，则管束中节点1的浓度即可得知，而管束中的节点2，3，4和5对应的浓度C₂，C₃，C₄和C₅分别为

$ \begin{array}{l} C_2^{n + 1} = \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_1^n + \left[ {1 - \frac{{2D\Delta t}}{{{{(\Delta x)}^2}}}} \right]C_2^n + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_3^n \end{array} $

(9)

$ \begin{array}{l} C_3^{n + 1} = \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_2^n + \left[ {1 - \frac{{2D\Delta t}}{{{{(\Delta x)}^2}}}} \right]C_3^n + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_4^n \end{array} $

(10)

$ \begin{array}{l} C_4^{n + 1} = \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_3^n + \left[ {1 - \frac{{2D\Delta t}}{{{{(\Delta x)}^2}}}} \right]C_4^n + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\frac{{D\Delta t}}{{{{(\Delta x)}^2}}} + \frac{{{v_m}\Delta t}}{{2\Delta x}}} \right]C_4^n \end{array} $

(11)

$ C_5^{n + 1} = \frac{{2D\Delta t}}{{{{(\Delta x)}^2}}}C_4^n + \left[ {1 - \frac{{2D\Delta t}}{{{{(\Delta x)}^2}}}} \right]C_5^n $

(12)

式中：C_i(i=1，2，3，4，5)为节点i的体积百分浓度，m³/m³；∆x为空间步长，m；∆t为时间步长，s；上标n为第n次迭代。

显式差分若要收敛和稳定必须对时间步长和空间步长加以限制。讷正^[23]的研究发现一维对流弥散方程显式差分必须满足式(13)和式(14)，本次研究所取参数均满足该条件。

$ {\Delta x < \frac{{2D}}{v}} $

(13)

$ {\Delta t < \frac{{{{(\Delta x)}^2}}}{{2D}}} $

(14)

2 网络模拟过程

基于混溶驱的基本理论及控制方程，按照以下7个步骤进行混溶驱的网络模拟。

(1) 构建多孔介质网络模型。按照多孔介质网络模型的建模方法^[22]构建模型大小为40×40×5的体中心(BCC)网络模型，管束半径服从对数均匀分布，不同的配位数通常通过随机断开管束的方式获得。不同参数的网络模型表征不同孔隙结构的多孔介质。

(2) 求解单相流动时网络模型压力场。在初始条件下，网络模型完全饱和被驱替流体，然后结合封闭边界条件和设定的驱替流量，求解模型的压力场分布，此时的压力场分布实际为被驱替流体单相流动时的压力分布。

(3) 确定计算时间步长∆T。根据被驱替流体单相流动的压力场分布计算流体通过单根管束的平均流体速度v_ij，Taylor-Aris弥散系数D和平均时间t_ij。

$ {{v_{ij}} = \frac{{{q_{ij}}}}{{{S_{ij}}}} = \frac{{{q_{ij}}}}{{\pi {\kern 1pt} r_{ij}^2}}} $

(15)

$ {{t_{ij}} = \frac{l}{{{v_{ij}}}}} $

(16)

式中：v_ij为通过单根管束ij的平均流体速度，m/s；S_ij为单根管束ij横截面积，m²；t_ij为通过单根管束ij的平均时间，s。

找出驱替流体通过单根管束的最短时间

$ {T_{{\rm{min}}}} = {\rm{min}}({t_{ij}}) $

(17)

式中：T_min为驱替流体通过单根管束的最短时间，s。

以最短时间的1/10为每次压力更新的时间步长(或计算时间步长)∆T

$ \Delta T = \frac{{{T_{{\rm{min}}}}}}{{10}} $

(18)

式中：∆T为计算时间步长，s。

计算时间步长∆T受模型非均质性、连通性以及驱替流量等因素的共同影响，变化范围大致为10^-4~10^-2 s。由于计算时间步长∆T较短，可认为在∆T内的驱替过程属于稳定驱替。差分时间步长∆t取为计算时间步长的1/100，即10^-6~10^-4 s，以保证显式差分格式收敛且稳定。

(4) 更新单管中的流体分布。不考虑弥散影响时，驱替流体在∆T内向前推进的距离为

$ \Delta {X_{ij}} = {v_{ij}}\Delta T $

(19)

驱替流体在单管中占据的总长度为所有∆X_ij的叠加，即X_ij =∑∆X_ij。若$\frac{l-X_{i j}}{v_{i j}} < \Delta T$，则须对计算时间步长∆T进行调整，∆T取为驱替前缘刚好到达管束节点的最短时间，即

$ \Delta T = {\rm{min}}\left( {\frac{{l - {X_{ij}}}}{{{v_{ij}}}}} \right) $

(20)

考虑弥散影响时，在所有与驱替流体接触的管束中，用式(9)~(12)的显式差分格式求解单根管束中各节点的浓度分布，结合流入节点的管束流量，利用加权平均的方法计算网络模型节点的浓度。

(5) 更新各管束传导率，为下一次压力求解作准备。不考虑弥散时，根据更新后驱替流体和被驱替流体所占管束长度更新管束传导率。考虑弥散时，单根管束为驱替流体和被驱替流体的混合带。由图 1可以发现，单根管束被剖分为4个网格，共使用了5个节点描述管束的浓度分布，每个网格的体积百分浓度可认为是网格两端点的平均值，因此，每一个网格的有效黏度为

$ {\mu _{ei}} = {\left[ {\frac{{{C_i} + {C_{i + 1}}}}{{2\mu _i^{0.25}}} + \frac{{1 - ({C_i} + {C_{i + 1}})}}{{2\mu _{i + 1}^{0.25}}}} \right]^{ - 4}} $

(21)

式中：μ_ei为网格i的有效黏度，mPa·s；μ_i，μ_{i + 1}分别为网格i，i + 1的流体黏度，mPa·s。

则式(3)也可被改写为

$ {q_{ij}} = \frac{{{{10}^9}\pi {\kern 1pt} r_{ij}^4({p_i} - {p_j})}}{{2l({\mu _{{\rm{e1}}}} + {\mu _{{\rm{e2}}}} + {\mu _{{\rm{e3}}}} + {\mu _{{\rm{e4}}}})}} $

(20)

(6) 根据修正后的有效黏度及传导率，利用共轭梯度法求解式(6)以得到定流量q所对应的压力场分布，随后更新各孔隙的平均流体速度和TaylorAris弥散系数D_ij。

(7) 重复步骤(4)，(5)和(6)实现流体的不断驱替，当流出物中驱替流体的相对浓度C/C₀达到0.999时结束模拟。

3 模拟结果

选取特定的网络模型参数(水力半径R = 4 μm、管束长度l = 200 μm、变异系数CV= 0.8、配位数z = 4.8)和流体运动参数(驱替流量q = 10-5 mL/s、被驱替流体与驱替流体的黏度比或流度比M = 10、分子扩散系数D_m = 10^-10 m²/s)作为混溶驱的标准条件，以此进行单因素分析。

3.1 不考虑和考虑弥散的混溶驱分析

弥散控制着驱替流体和被驱替流体的混合程度，不考虑弥散的混溶驱只存在对流作用，该过程与加入表面活性剂的油水非混溶驱类似，驱替过程没有考虑毛管压力的影响。

由图 2可看出：当被驱替流体与驱替流体的黏度比为1时，压差和压力场分布几乎不发生变化，当黏度比大于1时，压差随着注入孔隙体积(PV数)的增加而减小；在驱替刚开始时，考虑弥散的压差下降速度更慢，这有利于保持地层能量，扩大体积波及效率，提高采收率；当注入PV数达到0.8时，在不考虑弥散的情形中，驱替流体沿着贯穿的通道直接流出模型，被驱替流体很难再被驱替出来，压差趋于稳定，在考虑弥散的情形中，仍有部分被驱替流体被驱替出来，压差进一步减小；随着注入PV数的进一步增加，2种情形的压差都趋于一个定值，但是考虑弥散时对应的压差更小，主要是因为稳定情形类似于驱替流体的单相流动，考虑弥散的模型中驱替流体占据的体积相对更大，且驱替流体的黏度小于被驱替流体，导致最终压差偏小。

不考虑弥散时，驱替流体主要进入主流通道中^[24]，表现出黏性指进现象^[25][图 3(a)]。考虑弥散时，驱替流体以主流通道为基准向四周扩散，进入了更多的孔隙，其中包括部分流量较小的孔隙，从而导致沿主流通道流动的流体减少[图 3(b)]，虽然主流通道中的流体也会在弥散作用下发生向前的扩展，但是在该流量(v = 10^-5mL/s)的驱替下，其作用十分有限，并不能消除主流通道中流体减少带来的影响，因此驱替流体较为均匀的向前推进，压制了黏性指进现象。对比图 3(a)和图 3(b)可以发现，2种情形均是从同一根管束突破，说明多孔介质的无序性决定了流体流动的主流通道，而弥散作用只是将通道进行延展。

图 4为流出物中驱替流体相对浓度C/C₀随注入PV数的变化关系，该曲线类似于油水非混溶驱的含水率曲线。由图 4可以发现：不考虑弥散时，驱替流体在0.26 PV发生突破，对应的采出程度为26%，突破后，C/C₀迅速上升，达到0.7左右，然后再缓慢上升；考虑弥散时，驱替流体在0.36 PV发生突破，对应的采出程度为36%，比不考虑弥散对应的突破时间要晚，说明弥散作用有效压制了黏性指进现象，且C/C₀呈缓慢上升的形态，不存在骤变的情况，当注入PV数达到2时，C/C₀基本达到1。

3.2 考虑弥散的混溶驱影响因素分析 3.2.1 孔隙非均质性

图 5为标准条件下其他参数不变时，不同变异系数下驱替前缘突破时的二维图形。由图 5可以看出，C为0.05时，驱替流体均匀推进，且过渡带很短；随着CV增加，孔隙半径差异性增大，黏性指进现象越明显，形成的指进通道越多；随着CV进一步增加，模型中会出现难以被驱替的管束，驱替流体在此类管束附近发生绕流，呈现中心管束未被驱替而四周管束被驱替的情形，形成残余区域。

由图 6可以看出：当CV= 0.05时，前缘突破时间最晚，突破后C/C₀上升的速度最快，当注入PV数达到1.2时，C/C₀几乎不再发生变化，且等于1；随着CV的增加，前缘突破时间提前，而C/C₀达到1的时间延后，这大大增加了前缘突破后的生产时间；注入PV数为1所对应的C/C₀随CV的增加稍有增加。

由图 7可知：采出程度与CV满足线性关系，CV为0.05时对应的采出程度最高，为85%；随着CV增加，采出程度急剧下降，当CV达到1.05时，采出程度仅为12%。

3.2.2 孔隙连通性

孔隙连通性和孔隙非均质性都可以表征多孔介质微观孔隙结构的特征，二者对混溶驱替过程的影响是类似的，但孔隙非均质性对混溶驱的影响则更大。

总体来看，前缘突破时的采出程度与配位数(z)呈线性关系(图 8)。当模型完全连通(z = 8.0)时，前缘突破时间最晚，C/C₀到达1的时间最早，对应的采出程度最高，约为61%；随着z的减小，形成的指进通道增多，残余区域面积增大，前缘突破时间逐渐缩短，而C/C₀到达1的时间逐渐延长。

3.2.3 黏度比

被驱替流体和驱替流体的黏度比严重影响着混溶驱替过程。当黏度比(M)为2时，整个驱替过程基本呈活塞式驱替；随着M的增大，主流通道效应显著，整个驱替过程表现出黏性指进现象，且M越大，黏性指进现象越明显；当M达到100时，驱替流体几乎只沿一条通道流动，体积波及效率极低(图 9)。

随着M的增加，前缘突破时间逐渐缩短，而C/C₀到达1的时间逐渐延长(图 10)。当M增大到一定程度时(M ≥ 50)，前缘突破时间几乎不再发生变化，但前缘突破后C/C₀的上升速度和幅度受M影响，对于M为50的情形，C/C₀随着注入PV数的增加缓慢增加，对于M为100的情形，C/C₀快速上升到0.6左右，然后随着注入PV数的增加缓慢增加。

由图 11可看出：M越接近1，前缘突破时的采出程度越高；随着M的增加，采出程度急剧下降，其变化趋势基本服从乘幂规律；当M大于等于50时，前缘突破时的采出程度开始趋于平稳。由此可以发现，当M较大时，驱替流体会快速突破，使体积波及效率偏低，因此在CO₂驱油过程中常加入泡沫以增加流体在高渗通道中的渗流阻力或加入高分子量聚合物以增加CO₂黏度，以此防止CO₂快速突破。

3.2.4 驱替流量

驱替流量直接决定对流作用和弥散作用的强弱。当驱替流量较小时，对流作用相对较弱，弥散作用相对较强，驱替流体在弥散作用下可大范围向四周扩散，从而使驱替流体较为均匀地向前驱替，呈现活塞式驱替。随着驱替流量的增加，对流作用大大增强，而弥散作用稍有增强，导致驱替流体主要在主流通道流动，向四周扩散的范围有限，表现出黏性指进现象，且存在许多残余区域。当驱替流量增加到一定值(模拟条件下10^-4 mL/s)时，驱替过程几乎不再发生变化(图 12)。

在模拟条件下，驱替流量小于等于10^-6 mL/s时，驱替过程差异很小，前缘突破时间和C/C₀到达1的时间几乎一致，说明驱替流体受对流作用的影响十分有限，且流体在弥散作用下已经完成了所有的横向扩散，体积波及效率已经达到极致。随着驱替流量的增大，前缘突破时间提前，C/C₀到达1的时间延后。当驱替流量增大到一定程度时，对流作用将占主体地位，驱替流量对前缘突破时间的影响减弱，转为影响最终采收率(图 13)。

当注入流量小于等于10^-6 mL/s时，前缘突破时的采出程度基本保持在60%不变；当注入流量大于10^-6 mL/s时，采出程度随着注入流量的增加而减小，其变化趋势基本服从乘幂规律(图 14)。

4 结论

(1) 在体中心网络模型中通过对模型压力场和浓度场的耦合求解，提出了考虑弥散影响的混溶驱动态网络模拟方法。

(2) 多孔介质的无序性决定流体流动的主流通道，而弥散作用将主流通道进行延展。与不考虑弥散的混溶驱相比，考虑弥散的混溶驱的模型两端压差减小速度更慢，最终稳定时的压差更小，体积波及效率更大，驱替流体突破时间更晚，前缘突破时的采出程度更高，流出物中驱替流体的相对浓度C/C₀上升更为缓慢。

(3) 对于考虑弥散的混溶驱替过程，孔隙非均质性、孔隙连通性、黏度比和驱替流量等对混溶驱均具有较大影响，其中黏度比对混溶驱的影响相对较强，孔隙连通性对混溶驱的影响相对较弱。孔隙非均质性越强(或孔隙连通性越低或黏度比越大或驱替流量越大)，前缘突破时间越早，黏性指进现象越明显，体积波及效率越低。前缘突破时的采出程度与孔隙非均质性(或孔隙连通性)呈线性关系，与黏度比(或驱替流量)呈乘幂关系。