引言

在油田开发中后期，油藏压力降低，通常采用注水来维持地层压力并且提高采收率^[1-3]。然而，由于注入水和地层水不相容，会产生结垢沉淀，例如：CaSO$_4$·1/2H$_2$O，CaSO$_4$·2H$_2$O，CaSO$_4$，CaCO$_3$，BaSO$_4$，SrSO$_4$等。碳酸盐矿物形成后可以通过向地层中挤注酸来去除，但是硫酸盐矿物的防治仍旧是国内外的一大难题^[4-6]。

油田结垢预测模型主要有3类：(1) 定性结垢预测公式，以饱和指数为基础，主要有Rynar稳定指数法^[7]、D-S饱和指数法和Oddo-Tomson指数法^[8-9]；(2) 以热力学平衡原理为基础的热力学结垢模型，如pitzer离子相互作用模型^[10]、UNIQUAC通用似化学模型^[11]；(3) 基于动力学原理的结垢预测模型，Bedrikovetsky^[12-16]根据化学反应速率方程和物质的对流扩散方程建立了预测硫酸钡结垢的动力学模型。

以饱和指数为基础的定性结垢预测公式是经验公式，准确性不高；以热力学平衡原理建立的模型，没有考虑注水驱替时的动力学因素。本文以Bedrikovetsky的动力学模型为基础，在瞬态模型中考虑了扩散项，优化了求解模型参数的方法，并采用数值方法求解，得出驱替至不同PV时岩芯中结垢量的分布。

1 模型假设

(1) 注入水中的SO$_4^{2-}$和地层水中的Ba$^{2+}$的化学反应不可逆，为二级基元反应并遵循质量作用定律。

(2) 硫酸钡瞬时沉积，且生成的沉淀不随水流运移。

(3) 注入水与地层水不可压缩。

(4) 化学反应速率常数$K_{\rm{a}}$和扩散系数$D$均与流速$U$成比例，且所有物质的扩散系数相同。

(5) 岩芯的孔隙度ϕ和渗透率$K$分布均匀。

2 模型的建立及推导

影响硫酸钡结垢的因素主要有：对流、扩散和沉淀反应^[17]，因此在建立的控制方程中应该包括对流项、扩散项和反应项。注入水与地层水混合产生沉淀的流动控制方程应包括SO$_4^{2-}$、Ba$^{2+}$和BaSO$_4$的质量守恒方程^[18-20]。

$ \left\{ \begin{array}{l} \phi \dfrac{{\partial c_{{\rm{Ba}}^{2 + } } }}{{\partial t}} + U\dfrac{{\partial c_{{\rm{Ba}}^{2 + } } }}{{\partial x}} = D\dfrac{{\partial ^2 c_{{\rm{Ba}}^{2 + } } }}{{\partial x^2 }} - K_{\rm{a}} c_{{\rm{Ba}}^{2 + } } c_{{\rm{SO}}_4^{2 - } } \\[8pt] \phi \dfrac{{\partial c_{{\rm{SO}}_4^{2 - } } }}{{\partial t}} + U\dfrac{{\partial c_{{\rm{SO}}_4^{2 - } } }}{{\partial x}} = D\dfrac{{\partial ^2 c_{{\rm{SO}}_4^{2 - } } }}{{\partial x^2 }} - K_{\rm{a}} c_{{\rm{Ba}}^{2 + } } c_{{\rm{SO}}_4^{2 - } } \\[8pt] \phi \dfrac{{\rho _{_{\rm{BaSO_4}} } }}{{M_{_{{\rm{BaSO}}_4 } } }}\dfrac{{\partial \sigma }}{{\partial t}} = K_{\rm{a}} c_{{\rm{Ba}}^{2 + } } c_{{\rm{SO}}_4^{2 - } } \\ \end{array} \right. $

(1)

驱替开始时岩芯中已饱和富含Ba$^{2+}$的模拟地层水，因此初始条件可表示为

$ c_{{\rm{Ba}}^{2 + } } = c_{{\rm{Ba}}^{2 + } }^0 {\rm{, }}~~~c_{{\rm{SO}}_4^{2 - } } = 0 $

当富含SO$_4^{2-}$的模拟注入水瞬时注入岩芯内时，在岩芯入口端的SO$_4^{2-}$保持恒定，因此岩芯入口端的边界条件为

$ c_{{\rm{Ba}}^{2 + } } = 0{\rm{, }}~~~c_{{\rm{SO}}_4^{2 - } } = c_{{\rm{SO}}_4^{2 - } }^0 $

假设流体在岩芯出口流出，无倒流发生，则岩芯出口端的边界条件为$\dfrac{{\partial c_{{\rm{Ba}}^{2 + } } }}{{\partial x}} = \dfrac{{\partial c_{{\rm{SO}}_4^{2 - } } }}{{\partial x}} = 0$

为了消除模型中不同参数之间数量级的差异，需要首先定义无因次参数

$ c_{{\rm{Ba}}^{2+}}^{*}=\dfrac{{c_{{\rm{Ba}}^{{\rm{2 + }}} } }}{{c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 }}；{\rm{ }}c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^{*} = \dfrac{{c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} } }}{{c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^0 }}{\rm{ }}\\ {{S}} = \dfrac{{\rho_{_{\rm{BaSO_{\rm{4}}}} } }}{{M_{_{{\rm{BaSO}}_{\rm{4}} } } }}\dfrac{\sigma }{{c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 }}；{\rm{ }}c^0 = \dfrac{{c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 }}{{c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^0 }}；x_{\rm{d}} = \dfrac{x}{L}{\rm{ }}；$

$t_{\rm{d}}$；$ = \dfrac{{Ut}}{{\phi L}}{\rm{ }}$；$\varepsilon _{\rm{D}} = \dfrac{D}{{LU}}{\rm{ }}$；$\varepsilon _{\rm{k}} = \dfrac{{K_{\rm{a}} Lc_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} } }}{U} $对模型无因次化。

将定义的无因次参数代入式(1)，得到的无因次控制方程如式(2) 所示

$ \left\{ \begin{array}{l} \dfrac{{\partial c_{{\rm{Ba}}^{{\rm{2 + }}} }^* }}{{\partial t_{\rm{d}} }} + \dfrac{{\partial c_{{\rm{Ba}}^{{\rm{2 + }}} }^* }}{{\partial x_{\rm{d}} }} = \varepsilon _{\rm{D}} \dfrac{{\partial ^2 c_{{\rm{Ba}}^{{\rm{2 + }}} }^* }}{{\partial x_{\rm{d}} ^2 }} - \varepsilon _{\rm{k}} c_{{\rm{Ba}}^{{\rm{2 + }}} }^* c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* \\[8pt] \dfrac{{\partial c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* }}{{\partial t_{\rm{d}} }} + \dfrac{{\partial c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* }}{{\partial x_{\rm{d}} }} = \varepsilon _{\rm{D}} \dfrac{{\partial ^2 c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* }}{{\partial x_{\rm{d}} ^2 }} - \varepsilon _{\rm{k}} c^0 c_{{\rm{Ba}}^{{\rm{2 + }}} }^* c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* \\[8pt] \dfrac{{\partial S}}{{\partial t_{\rm{d}} }} = \varepsilon _{\rm{k}} c_{{\rm{Ba}}^{{\rm{2 + }}} }^* c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* \\[8pt] t_{\rm{d}} = 0:{\rm{ }}c_{{\rm{Ba}}^{{\rm{2 + }}} }^* = 1, {\rm{ }}c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* = 0 \\[8pt] t_{\rm{d}} > 0, {\rm{ }}x_{\rm{d}} = 0:{\rm{ }}c_{{\rm{Ba}}^{{\rm{2 + }}} }^* = 0, {\rm{ }}c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* = 1 \\[8pt] t_{\rm{d}} > 0, x_{\rm{d}} = L:\dfrac{{\partial c_{{\rm{Ba}}^{{\rm{2 + }}} }^* }}{{\partial x_{\rm{d}} }} = \dfrac{{\partial c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* }}{{\partial x_{\rm{d}} }} = 0 \\ \end{array} \right. $

(2)

通过以上的控制方程、初始条件和边界条件，可以求出SO$_4^{2-}$、Ba$^{2+}$的浓度在岩芯中分布的数值解，进而得出岩芯中不同驱替PV时岩芯中硫酸钡沉淀量的变化。

3 硫酸钡结垢量预测

选取2块长庆油田的岩芯，岩芯物性见表 1。

实验使用的注入水和地层水分别为洛河水和H51-20井地层水，注入水中含有较多的SO$_4^{2-}$(成垢阴离子)，矿化度为4 737.01 mg/L，水型为Na$_2$SO$_4$型；地层水含有较多的Ba$^{2+}$(成垢阳离子)，矿化度为4 674.03 mg/L，水型为CaCl$_2$型。为了避免发生其他成垢离子之间的化学反应，在保证矿化度不变的前提下去除了注入水和地层水中的其他成垢离子，水样的离子成分如表 2所示。

由表 2可知，地层水中的Ba$^{2+}$与注入水中的SO$_4^{2-}$的初始浓度比$c^0 = 0.550 25$。

岩芯驱替实验的设备和流程如图 1所示。

3.1 扩散系数$D$的求取

设计岩芯驱替实验，用去除SO$_4^{2-}$的模拟注入水驱替饱和模拟地层水的岩芯，测出驱替不同PV时岩芯出口端Ba$^{2+}$的浓度来拟合扩散系数$D$。该过程可以用对流扩散方程来描述，见方程(3)

$ \left\{ \begin{array}{l} \phi \dfrac{{\partial c_{{\rm{Ba}}^{{\rm{2 + }}} } }}{{\partial t}} + U\dfrac{{\partial c_{{\rm{Ba}}^{{\rm{2 + }}} } }}{{\partial x}} = D\dfrac{{\partial ^2 c_{{\rm{Ba}}^{{\rm{2 + }}} } }}{{\partial x^2 }} \\ t = 0:{\rm{ }}c_{{\rm{Ba}}^{{\rm{2 + }}} } = \left\{ \begin{array}{l} 0, {\rm{ }}x \leqslant 0 \\ c_{{\rm{Ba}}^{{\rm{2 + }}} }^0, {\rm{ }}x > 0 \\ \end{array} \right. \\ c_{{\rm{Ba}}^{{\rm{2 + }}} } (x \to - \infty, t) = 0 \\ c_{{\rm{Ba}}^{{\rm{2 + }}} } (x \to + \infty, t) = c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 \\ \end{array} \right. $

(3)

微分方程(3) 的解析解如式(4) 所示

$ c_{{\rm{Ba}}^{{\rm{2 + }}} } (x, t) = c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 - \dfrac{{c_{{\rm{Ba}}^{{\rm{2 + }}} }^0 }}{2}{\rm{erfc}}\left( {\dfrac{{x - Ut/\phi }}{{2\sqrt {Dt/\phi } }}} \right) $

(4)

式(4) 的无因次形式可以表示为

$ c_{{\rm{Ba}}^{{\rm{2 + }}} }^* (x_{\rm{d}}, t_{\rm{d}} ) = 1 - \dfrac{1}{2}{\rm{erfc}}\left( {\dfrac{{x_{\rm{d}} - t_{\rm{d}} }}{{2\sqrt {\varepsilon _{\rm{D}} t_{\rm{d}} } }}} \right) $

(5)

式(5) 在$x_{\rm{d}} = 0$时不连续的物理初始条件可以通过连续函数$c_{{\rm{Ba}}^{{\rm{2 + }}} } (x, t)$在$t_{\rm{d}} \to 0$时的极限来模拟。如图 2所示，当$t_{\rm{d}}$趋近于0时，$x_{\rm{d}}$=0处的无因次钡离子浓度无限趋近于1，即此时岩芯入口端的钡离子浓度为初始钡离子浓度，即可认为岩芯入口端的初始条件连续。

因此，驱替实验中岩芯出口端的无因次Ba$^{2+}$的浓度为

$ c_{{\rm{Ba}}^{{\rm{2 + }}} }^* (x_{\rm{d}}, t_{\rm{d}} ) = 1 - \dfrac{1}{2}{\rm{erfc}}\left( {\dfrac{{x_{\rm{d}} - t_{\rm{d}} }}{{2\sqrt {\varepsilon _{\rm{D}} t_{\rm{d}} } }}} \right) $

(6)

设定驱替流速为0.5 mL/min，在驱替至不同PV时记录岩芯出口端的Ba$^{2+}$浓度，见表 3。实验中注入1 PV的时间$t_{{\rm{1PV}}} = \dfrac{{V_{\rm{P}} }}{{AU}} = \dfrac{{\phi L}}{U}$，无因次时间$t_{\rm{d}} = \dfrac{{Ut_{{\rm{1PV}}} }}{{\phi L}} = 1$，因此，注入1 PV的时间恒等于无因次时间$t_{\rm{d}}$。

由表 3中的岩芯出口端无因次Ba$^{2+}$浓度结合式(6) 可以分别拟合出两块岩芯的无因次扩散系数$\varepsilon _{\rm{D}}$，见图 3和图 4。

由图 3和图 4可知，当不发生结垢反应时，随着无因次驱替时间$t_{\rm{d}}$的增加，岩芯出口端的Ba$^{2+}$浓度逐渐减小；驱替2 PV后，1号与2号岩芯出口端的Ba$^{2+}$浓度已分别降低至初始Ba$^{2+}$浓度的11.15%和8.88%；经式(6) 拟合可得，1号岩芯的无因次扩散系数为0.168 3，2号岩芯的无因次扩散系数为0.137 6。

由式(2) 可得1号岩芯和2号岩芯的扩散系数，见表 4。

3.2 反应速率常数$K_{\rm{a}}$的求取

求取反应速率常数的一种方法是直接对式(2) 求解析解，得出稳态条件下Ba$^{2+}$和SO$_4^{2-}$在岩芯空间的分布，再通过实验拟合出无因次反应速率常数，但是这种方法需要在假设$c^0\ll1$的条件下进行。因此，该方法在应用时具有较大的局限性。本文改进了反应速率常数$K_{\rm{a}}$的拟合方法，直接求式(2) 在区域$0 \leqslant x_{\rm{d}} \leqslant 1, 0 \leqslant t_{\rm{d}} \leqslant t_{{\rm{df}}}$上的数值解，并结合3.1中得出的无因次扩散系数，对岩芯出口端的成垢离子浓度进行预测，再利用驱替实验测出的岩芯出口端离子浓度，反演出无因次反应速率常数$\varepsilon _{\rm{k}}$，提高了模型的适用性。

在求取扩散系数时驱替实验中的注入水中不含成垢阴离子，因此岩芯中无结垢发生。将岩芯1和岩芯2再次饱和模拟地层水，利用含有SO$_4^{2-}$的模拟注入水驱替饱和模拟地层水的岩芯，模拟注入水中的SO$_4^{2-}$与模拟地层水中的Ba$^{2+}$发生化学反应，设定流速为0.5 mL/min，在岩芯出口端测定驱替不同PV下的Ba$^{2+}$浓度，实验数据见表 5。

定义$X_j = c_{{\rm{Ca}}^{{\rm{2 + }}} }^* (x_{{\rm{d}}j}, t_{{\rm{d}}n} )$，$Y_j = c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^* (x_{{\rm{d}}j}, t_{{\rm{d}}n} )$，因此式(2) 可以表述为代数方程式的形式

$ \left\{ \begin{array}{l} \dfrac{{{\rm{d}}X_j }}{{{\rm{d}}t_{\rm{d}} }} = \varepsilon _{\rm{D}} \dfrac{{X_{j - 1} - 2X_j + X_{j + 1} }}{{(\Delta x_{\rm{d}} )^2 }} - \dfrac{{X_{j + 1} - X_j }}{{\Delta x_{\rm{d}} }} - \varepsilon _k X_j Y_j {\rm{ }}~~~~~j{\rm{ = 2, 3, }}...N \\[10pt] \dfrac{{{\rm{d}}Y_j }}{{{\rm{d}}t_{\rm{d}} }} = \varepsilon _D \dfrac{{Y_{j - 1} - 2Y_j + Y_{j + 1} }}{{(\Delta x_{\rm{d}} )^2 }} - \dfrac{{Y_{j + 1} - Y_j }}{{\Delta x_{\rm{d}} }} - \varepsilon _k c^0 Y_j X_j {\rm{ }}~~~~~j{\rm{ = 2, 3, }}...N \\ \dfrac{{\partial S}}{{\partial t_{\rm{d}} }} = \varepsilon _k X_j Y_j \\ t_{\rm{d}} = 0, {\rm{ }}X_j = 1, Y_j = 0 \\ j = 1, {\rm{ }}X_1 = 0, Y_1 = 1 \\ j = N + 1, {\rm{ }}X_{N + 1} = X_N, Y_{N + 1} = Y_N \\ \end{array} \right. $

(7)

最后，这一初值问题的方程组在$0 \leqslant t_{\rm{d}} \leqslant t_{{\rm{df}}}$的范围内可以运用第四和第五阶Runge-Kutta方程求解。求解结果如图 5和图 6所示。

由图 5和图 6可知，当含有SO$_4^{2-}$的注入水驱替饱和地层水的岩芯后，由于发生了结垢反应，岩芯出口端Ba$^{2+}$浓度迅速下降，SO$_4^{2-}$浓度迅速增大，呈现出“S”型曲线的趋势；当驱替至1 PV时，1号和2号岩芯出口端的Ba$^{2+}$浓度分别为初始Ba$^{2+}$浓度的1.14%和2.51%，即此时岩芯内的Ba$^{2+}$已经几乎被消耗殆尽；经式(7) 拟合可得1号岩芯的无因次反应速率常数为3.372 9，2号岩芯的无因次反应速率常数为1.867 5。

把无因次反应速率常数代入定义的无因次量中，即可求出岩芯1和岩芯2的反应速率常数，见表 6。

3.3 硫酸钡结垢量预测

在得出模型的无因次反应速率常数和无因次扩散系数之后，对式(1) 在$0 \leqslant x \leqslant L$，$0 \leqslant t \leqslant t_{\rm{f}} $上进行差分，便可得出驱替至不同PV时岩芯内SO$_4^{2-}$和Ba$^{2+}$的浓度变化，进而可得出驱替至不同PV时岩芯内硫酸钡结垢量的变化，结果见图 7和图 8。

由图 7和图 8可知，当含有SO$_4^{2-}$的模拟注入水注入岩芯后，两块岩芯内均发生了不同程度的结垢现象。随着模拟注入水注入体积的增加，岩芯内硫酸钡结垢量均逐渐增大；最大结垢量的位置向岩芯中部转移；当注入孔隙体积倍数小于1 PV时，岩芯内硫酸钡结垢量迅速增大，当注入孔隙体积倍数大于1 PV时，硫酸钡结垢量变化不大；在沿岩芯长度方向上，硫酸钡结垢量呈现先增大后减小的趋势。

4 结论

(1) 通过测定驱替实验中岩芯出口端Ba$^{2+}$的浓度并对模型求解，反演出了模型中的两个重要参数——扩散系数和反应速率常数。

(2) 对含扩散项的硫酸钡结垢动力学模型有限差分并采用第四阶和第五阶Runge-Kutta方程求解，可以避免求解解析解时岩芯入口端Ba$^{2+}$与SO$_4^{2-}$浓度比的假设。

(3) 岩芯出口端的Ba$^{2+}$浓度随着驱替时间的增加而降低，硫酸根离子浓度随着驱替时间的增大而升高。

(4) 岩芯中的硫酸钡结垢量随模拟注入水注入量的增大而增大。

(5) 随着模拟注入水注入体积的增加，最大结垢量的位置逐步向岩芯中部转移。

(6) 在沿岩芯长度方向上，结垢量呈现出先增加后减小的趋势。

符号说明

ϕ—岩芯孔隙度，无因次；

$c_{{\rm{Ba}}^{2 + } }$—钡离子浓度，mg/L；

$t$—时间，s；

$U$—流速，m/s；

$D$—扩散系数，m$^2$/s；

$K_{\rm{a}}$—反应速率常数，mol/(L·s$^{-1}$)；

$c_{{\rm{SO}}_4^{2 - } }$—硫酸根离子浓度，mg/L；

$\rho$—密度，kg/m$^3$；

$M$—摩尔质量，kg/mol；

$\sigma$—硫酸钡结垢量，mg/L；

$c_{{\rm{Ba}}^{2+}}^{*}$—无因次钡离子浓度；

$c_{{\rm{Ba}}^{{\rm{2 + }}} }^0$—初始钡离子浓度，mg/L；

$c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^{*}$—无因次硫酸根离子浓度；

$c_{{\rm{SO}}_{\rm{4}}^{{\rm{2 - }}} }^0$—初始硫酸根离子浓度，mg/L；

$S$—无因次硫酸钡结垢量；

$c^0$—初始钡离子浓度与硫酸根离子浓度比；

$x_{\rm{d}}$—无因次岩芯长度；

$t_{\rm{d}}$—无因次时间；

$L$—岩芯长度，cm；

$\varepsilon _{\rm{D}}$—无因次扩散系数；

$\varepsilon_{\rm{k}}$—无因次反应速率常数；

$p$—压力，Pa；

$t_{\rm{df}}$—无因次时间右边界，也可以理解为限定的无因次时间，无因次；

$\mu$—黏度，mPa·s。