·726· 工程科学学报，第41卷，第6期 X.+X。=1 求解步骤：(1)计算不同矿化度下黏土扩散双 (3)初始条件. 电层的电势、微粒的临界滞留体积分数及迁移体 Pl=0=PiniS=0=SveCI=0=Cimi 积分数：(2)计算黏土微粒运移对水相渗透率的影 (4)边界条件. 响，绘制矿化度与水相渗透率变化倍数曲线：（3) 9,lx=0=9.9,1x=L=-9q.=9。+9。+9e MPES方法求解n时刻的压力、含水饱和度分 式中，x为距离，m;t为时间，s;P为压力，MPa:S、S。 布B、显式求解地层水矿化度分布：(4)根据步骤 分别为含水饱和度和含油饱和度；u。以.分别为油 (2)的结果，插值计算n时刻各网格渗透率损伤倍 和水的黏度，mPa·s;X。、X.分别为盐、水组分在水相 数，更新各网格渗透率；(5)判断是否小于模拟时 中所占的质量分数；S为束缚水饱和度；Cm为油藏 间，是则返回步骤(3)计算n+1时刻的各参数，否 初始地层水矿化度，mgL;9.,q。和q。分别为水、 则循环结束 油、盐组分的体积流量，m3·d1.则以上油、水、盐相 3结果与分析 方程及辅助方程、初始条件和边界条件构成了耦合 了微粒运移机理和储层伤害机理的完整低矿化度水3.1离子矿化度与价型对电势的影响 驱数学模型 以NaCl溶液为例，利用式(4)、式(5)计算不同 渗流模型中，水相中盐组分所占的体积比重很 矿化度(58.4、584和5840mgL-)下扩散双电层厚 小，所以水相方程求解时X≈1，联立油、水的连续 度和Zeta电势与表面电势的比值变化规律曲线，如 性方程，结合辅助方程将含水饱和度消掉，得到只含 图5所示.同理计算矿化度为292mgL-1时离子价 有压力的偏微分方程： 型分别为l、2、3的扩散双电层厚度和Zeta电势与 云a+两u)+是（低g) a KK ap 表面电势的比值变化规律曲线，如图6所示 k。成+g.=0 1.0 (28) 0.9 ◆C=58.4mgL ★一C=584mgL 模型离散为： 一C=5840mg·L Pi-PI P-PRI 0.7 入 △x -入-1 △x △x -+9.=0(29) 结合内外边界条件得到方程组，采用追赶法进 03 行求解，即可得到压力分布.然后将求得的压力带 入水相的连续性方程中，可得差分方程： 02 0.1 Pi-P --入-1 P-P 0 △x △x 10 20 30 40 50 扩散双电层厚度mm △x -+9m= 6- 图5不同矿化度下的Zta电势与表面电势比值分布曲线 (30) Fig.5 Curves of Zeta potential to surface potential ratio at a different △L ion concentrations A= KiKo (31) 4.(1+BG) 由图5~6可以看出，双电层性质对离子矿化度 式中，入”为n时刻网格i处水相的流动系数， 和价型较为敏感.随着离子矿化度和价型数的降 mD·(mPas)-:△t为差分过程中的时间步长，s;下 低，双电层膨胀变厚，（电势升高，（电势的升高破坏 标i为网格编号；上标n为时间编号. 了岩石表面黏土颗粒的受力平衡，黏土微粒失稳并 含水饱和度可以通过(30)式进行显式求解.假 发生运移.这是低矿化度水驱微粒运移机理提高采 设只考虑注入水与地层水之间的反应造成的离子矿 收率的根本原因. 化度的变化，盐组分的连续性方程可以简化为： 3.2基于DLVO理论的势能和静电力计算 CiniSwe+CiSi+Cini (Smi-S)= 取黏土微粒半径0.5μm,以NaCl溶液为例，计 Ca+1.S (32) 算离子价型为z=1时不同离子矿化度(58.4、584 式中：Ca为初始注入水矿化度，mgL;C:为n时 和5840mgL-)下的势能曲线，如图7所示：同理计 刻第i个网格的矿化度，mgL1. 算矿化度为292mgL-'时离子价型分别为1、2、3的工程科学学报，第 41 卷，第 6 期 Xw + Xc = 1 ( 3) 初始条件． P | t = 0 = Pini ; Sw | t = 0 = Swc ; C | t = 0 = Cini ( 4) 边界条件． qv | x = 0 = qv ; qv | x = L = － qv ; qv = qw + qo + qc 式中，x 为距离，m; t 为时间，s; P 为压力，MPa; Sw、So 分别为含水饱和度和含油饱和度; μo、μw 分别为油 和水的黏度，mPa·s; Xc、Xw 分别为盐、水组分在水相 中所占的质量分数; Swc为束缚水饱和度; Cini为油藏 初始地层水矿化度，mg·L － 1 ; qw，qo 和 qc 分别为水、 油、盐组分的体积流量，m3 ·d － 1 ． 则以上油、水、盐相 方程及辅助方程、初始条件和边界条件构成了耦合 了微粒运移机理和储层伤害机理的完整低矿化度水 驱数学模型． 渗流模型中，水相中盐组分所占的体积比重很 小，所以水相方程求解时 Xw≈1，联立油、水的连续 性方程，结合辅助方程将含水饱和度消掉，得到只含 有压力的偏微分方程:   ( x KKrw ( 1 + βσstr) μw P  ) x +   ( x KKro μo P  ) x + qv = 0 ( 28) 模型离散为: λi ·Pn + 1 i + 1 － Pn + 1 i Δx － λi － 1·Pn + 1 i － Pn + 1 i － 1 Δx Δx + qv = 0 ( 29) 结合内外边界条件得到方程组，采用追赶法进 行求解，即可得到压力分布． 然后将求得的压力带 入水相的连续性方程中，可得差分方程: λn w·i Pn + 1 i + 1 － Pn + 1 i Δx － λn wi － 1·Pn + 1 i － Pn + 1 i － 1 Δx Δx + qwvi =  Sn + 1 wi － Sn wi Δt ( 30) λn wi = Kn i Krw μw ( 1 + βσstr) ( 31) 式中，λn wi 为 n 时 刻 网 格 i 处水相的流动系数， mD·( mPa·s) － 1 ; Δt 为差分过程中的时间步长，s; 下 标 i 为网格编号; 上标 n 为时间编号． 含水饱和度可以通过( 30) 式进行显式求解． 假 设只考虑注入水与地层水之间的反应造成的离子矿 化度的变化，盐组分的连续性方程可以简化为: Cini·Swc + Cn i Sn wi + Cinj ( Sn + 1 wi － Sn wi ) = Cn + 1 i ·Sn + 1 wi ( 32) 式中: Cinj为初始注入水矿化度，mg·L － 1 ; Cn i 为 n 时 刻第 i 个网格的矿化度，mg·L － 1 ． 求解步骤: ( 1) 计算不同矿化度下黏土扩散双 电层的 ζ 电势、微粒的临界滞留体积分数及迁移体 积分数; ( 2) 计算黏土微粒运移对水相渗透率的影 响，绘制矿化度与水相渗透率变化倍数曲线; ( 3) IMPES 方 法 求 解 n 时 刻 的 压 力、含 水 饱 和 度 分 布［31］、显式求解地层水矿化度分布; ( 4) 根据步骤 ( 2) 的结果，插值计算 n 时刻各网格渗透率损伤倍 数，更新各网格渗透率; ( 5) 判断是否小于模拟时 间，是则返回步骤( 3) 计算 n + 1 时刻的各参数，否 则循环结束． 3 结果与分析 3. 1 离子矿化度与价型对 ζ 电势的影响 以 NaCl 溶液为例，利用式( 4) 、式( 5) 计算不同 矿化度( 58. 4、584 和 5840 mg·L － 1 ) 下扩散双电层厚 度和 Zeta 电势与表面电势的比值变化规律曲线，如 图 5 所示． 同理计算矿化度为 292 mg·L － 1时离子价 型分别为 1、2、3 的扩散双电层厚度和 Zeta 电势与 表面电势的比值变化规律曲线，如图 6 所示． 图 5 不同矿化度下的 Zeta 电势与表面电势比值分布曲线 Fig． 5 Curves of Zeta potential to surface potential ratio at a different ion concentrations 由图 5 ～ 6 可以看出，双电层性质对离子矿化度 和价型较为敏感． 随着离子矿化度和价型数的降 低，双电层膨胀变厚，ζ 电势升高，ζ 电势的升高破坏 了岩石表面黏土颗粒的受力平衡，黏土微粒失稳并 发生运移． 这是低矿化度水驱微粒运移机理提高采 收率的根本原因． 3. 2 基于 DLVO 理论的势能和静电力计算 取黏土微粒半径 0. 5 μm，以 NaCl 溶液为例，计 算离子价型为 z = 1 时不同离子矿化度( 58. 4、584 和 5840 mg·L － 1 ) 下的势能曲线，如图7 所示; 同理计 算矿化度为 292 mg·L － 1时离子价型分别为 1、2、3 的 · 627 ·