交联聚合物驱是在聚合物驱基础上发展起来的一项提高原油采收率新技术^[1-9]。油藏数值模拟是交联聚合物驱研究的一个重要手段^[10-14]，模拟结果可为交联聚合物驱的驱油机理分析、注采条件优化、矿场实施效果预测和评价提供决策依据。在数值模拟过程中，超过80%的计算时间都用于求解渗流微分方程差分离散后得到的大型稀疏线性代数方程组，其性能好坏决定了油藏模拟的速度和质量^[15-16]。其中，渗流微分方程差分格式的选择直接影响求解的稳定性和计算效率。为此，本文提出交联聚合物驱渗流微分方程采用Douglas差分格式，利用分步降维的思路并结合追赶法进行求解，能够大幅度提高模型的运算速度，同时保证了运算精度及稳定性。

1 交联聚合物驱数学模型

交联聚合物驱是一个复杂的化学驱油过程，为了便于研究并尽可能多地考虑交联聚合物驱油机理，作如下假设:(1) 油藏中存在油、水两相渗流，油相中仅包含油组分，其他组分在水相中；(2) 流体渗流符合达西定律，且油、水相之间无质量交换；(3) 油藏岩石及流体微可压缩；(4) 考虑油、水相之间毛管压力和重力的影响；(5) 聚合物、交联剂等物质的加入不影响水相体积和密度的变化，不影响油相渗流和油相性质的变化；(6) 组分物质的扩散遵循Fick扩散定律，吸附遵循Langmuir等温吸附理论；(7) 化学反应只考虑聚合物和交联剂之间的交联反应；(8) 整个过程为等温渗流，因化学反应引起的温度变化忽略不计。在此基础上，推导得到交联聚合物驱的数学模型。

油组分的基本渗流微分方程

$\nabla\!\cdot\!\left[\!{{F_{{\rm{Dc}}}}\!\cdot \dfrac{{{K_{\rm{o}}}}}{{{B_{\rm{o}}}{\mu _{\rm{o}}}}}\nabla\!\left(\! {{p_{\rm{o}}}\!-\!{F_{\rm{g}}} \cdot {\rho _{\rm{o}}}{\rm{g}}D} \right)} \right]\!+\! {q_{{\rm{ov}}}}\!=\!\dfrac{ \partial }{{ \partial t}}\left(\! {\dfrac{{\phi\!{S_{\rm{o}}}}}{{{B_{\rm{o}}}}}}\!\right) $

(1)

水组分的基本渗流微分方程

$\nabla\!\cdot \left[\!{{F_{{\rm{Dc}}}}\!\cdot\! \dfrac{{{K_{\rm{w}}}}}{{{R_{{\rm{Kw}}}}{B_{\rm{w}}}{\mu _{\rm{w}}}}}\nabla \left(\!{{p_{\rm{w}}}\!-\!{F_{\rm{g}}}\!\cdot {\rho\!_{\rm{w}}}{\rm{g}}D} \right)} \right]\!+\!{q_{{\rm{wv}}}}\!=\hspace{2em}\!\dfrac{ \partial }{{ \partial t}}\left(\! {\dfrac{{\phi\!{S_{\rm{w}}}}}{{{B_{\rm{w}}}}}}\!\right)$

(2)

聚合物、交联聚合物和交联剂的基本渗流微分方程

$\begin{align} & \nabla \cdot \left[ {{F}_{\text{Dc}}}\cdot \frac{{{K}_{\text{w}}}}{{{R}_{\text{Kw}}}{{\mu }_{\text{w}}}}\nabla \left( {{p}_{\text{w}}}-{{F}_{\text{g}}}\cdot {{\rho }_{\text{w}}}\text{g}D \right){{C}_{l\text{w}}}+\phi {{S}_{\text{w}}}{{\beta }_{l}}\nabla {{C}_{l\text{w}}} \right] \\ & +{{B}_{\text{w}}}{{q}_{\text{wv}}}{{C}_{l\text{w}}}+\phi {{S}_{\text{w}}}{{R}_{l}}=\frac{\partial }{\partial t}\left( \phi {{E}_{l}}{{S}_{\text{w}}}{{C}_{l\text{w}}}+\phi {{C}_{l\text{ad}}} \right) \\ \end{align}$

(3)

毛管压力方程

${p_{{\rm{cow}}}} = {p_{\rm{o}}} - {p_{\rm{w}}} = {p_{\rm{c}}}\left( {{S_{\rm{w}}},{\sigma _{{\rm{wo}}}}} \right)$

(4)

饱和度归一化方程

${S_{\rm{o}}} + {S_{\rm{w}}} = 1$

(5)

2 物化机理表征

为了简单、准确地表征交联聚合物驱过程中的各种物理化学现象，对以下主要物化参数进行描述。

2.1 交联聚合物溶液黏度

计算交联聚合物溶液黏度随浓度变化的方法较多，大多是一些经验关系式。考虑到一种交联聚合物体系对应一组特定的黏浓关系曲线。本文采用液相黏度的非线性混合函数插值关系式来确定水相黏度。

$\ln{\mu_{\rm{w}}}=f({x_{\rm{P}}})\ln{\mu_{\rm{P}}}+ f({x_{\rm{G}}})\ln {\mu _{\rm{G}}}+\hspace{2em}\dfrac{{1-f({x_{\rm{P}}})- f({x_{\rm{G}}})}}{{1-{x_{\rm{P}}}-{x_{\rm{G}}}}}{x_{\rm{w}}}\ln\mu _{\rm{w}}^{'}$

(6)

2.2 吸附滞留

利用Langmuir吸附等温式计算组分l的吸附浓度大小${C_{l{\rm{ad}}}}$

${C_{l{\rm{ad}}}} = \dfrac{{a{C_{l{\rm{w}}}}}}{{1 + b{C_{l{\rm{w}}}}}}$

(7)

2.3 渗透率降低系数

由于聚合物、交联聚合物组分的吸附滞留，使得岩石孔隙介质的局部渗透率降低，定义水相渗透率降低系数来描述这一物理化学现象，计算关系式为

${R_{l{\rm{Kw}}}} = 1 + \dfrac{{\left( {{F_{l{\rm{rr}}}} - 1} \right){C_{l{\rm{ad}}}}}}{{{C_{l{\rm{ad}}\max }}}}$

(8)

${K_{\rm{w}}} = {K_{\rm{w}}^0}/\prod\limits_{l = {\rm{P}}，{\rm{G}}} {{R_{l{\rm{Kw}}}}} $

(9)

2.4 可及孔隙体积

聚合物和交联聚合物组分是高分子化合物，流经岩石多孔介质时不能遍及全部孔隙体积，将组分物质能到达的孔隙体积定义为可及孔隙体积。假定组分物质可及孔隙体积大小不变，岩石孔隙度为$\phi$，聚合物和交联聚合物组分能够到达部分的孔隙度分别为${\phi _{\rm{P}}}$和${\phi _{\rm{G}}}$ ，则可及孔隙体积倍数$E_{l}$可定义为

$E_l = \dfrac{{\phi _l }}{\phi },l = {\rm{G}}，{\rm{P}}$

(10)

2.5 扩散系数

聚合物、交联聚合物或交联剂组分的扩散是通过定义该组分在水相中的扩散系数来体现的，其扩散通量与扩散系数间的关系可通过Fick扩散定律进行描述

${F_{lJ}} = - {D_{lJ}}\dfrac{{ \partial {C_{l{\rm{w}}}}}}{{ \partial J}}\quad \left\{ {\begin{array}{*{20}{l}} {l = {\rm{P}}，{\rm{G}}，{\rm{Cr}}}\\ {J = x，y，z} \end{array}} \right.$

(11)

2.6 交联反应动力学

目前，聚合物交联反应体系主要包括过渡金属交联体系、醛类交联体系和复合交联体系3大类^[17]。其中，过渡金属交联体系中的聚合物/氯化铬模型最常见，该模型的反应方程式简记为

${\rm{polymer}} + n{\rm{Cr}}^{{\rm{3}} + } + m{\rm{H}}_2 {\rm{O}} \to {\rm{clpolymer}} + 3n{\rm{H}}^ +$

(12)

交联反应动力学模型为

$\dfrac{{{\rm{d}}{C_{{\rm{Cr}}}}}}{{{\rm{d}}t}} = - n \cdot \alpha \cdot {C_{{\rm{Cr}}}}^{\rm{2}} \cdot {C_{\rm{P}}}^{\rm{2}}$

(13)

$ - \dfrac{{{\rm{d}}{C_{\rm{P}}}}}{{{\rm{d}}t}} = \dfrac{{{\rm{d}}{C_{\rm{G}}}}}{{{\rm{d}}t}} = \alpha \cdot {C_{{\rm{Cr}}}}^{\rm{2}} \cdot {C_{\rm{P}}}^{\rm{2}}$

(14)

3 模型求解

模型求解采用隐式求解压力、显式求解饱和度与组分浓度的方法。

对于二维或三维渗流问题的压力方程，通过常规的差分格式 (如Crank-twodash Nicolson差分格式)离散后得到的是一个五对角或七对角的线性方程组，随着网格数增加，运算量将呈级数倍增长，即使选用各种迭代算法，也很难降低计算量。为提高运算效率和计算精度，提出采用计算量较小而又无条件稳定的Douglas差分格式对交联聚合物驱数学模型进行求解。

1955年，Peaceman D W和Rachford H H提出了计算量小、具有二阶误差精度且绝对稳定的P-twodash R差分格式^[18]，取得了较好的计算效果并受到人们的普遍关注，但它的缺点是不能推广到三维情形(直接推广到三维是绝对不稳定的)。1962年，Douglas J提出的Douglas差分格式^[19]克服了P-R格式的缺点，与P-R差分格式一样，每一时间步的计算可分解成若干个一维隐格式，只要求解具有三对角系数矩阵的方程组，从而可用追赶法求解。

假设三维抛物型方程形式如下

$ \dfrac{{{ \partial ^2}p}}{{ \partial {x^2}}} + \dfrac{{{ \partial ^2}p}}{{ \partial {y^2}}} + \dfrac{{{ \partial ^2}p}}{{ \partial {z^2}}} = \dfrac{{ \partial p}}{{ \partial t}}$

(15)

对应于式(15)的Douglas差分格式为

$ \left\{ {\begin{array}{*{20}{l}} {\dfrac{{p_{i,j,k}^{n+\frac{1}{3}}\!-\!p_{i,j,k}^n}}{{\Delta t}}\!=\! \dfrac{1}{{{h^2}}}\left[\!{\delta _x^2\left(\!{\dfrac{{p_{i,j,k}^{n+ \frac{1}{3}}\!+\!p_{i,j,k}^n}}{2}}\!\right)\!+\!\left(\!{\delta _y^2\! +\!\delta\!_z^2}\!\right)p_{i,j,k}^n}\!\right]}\\ {\dfrac{{p_{i,j,k}^{n + \frac{2}{3}} - p_{i,j,k}^{n + \frac{1}{3}}}}{{\Delta t}} = \dfrac{1}{{{h^2}}}\delta _y^2\left( {\dfrac{{p_{i,j,k}^{n + \frac{2}{3}} - p_{i,j,k}^n}}{2}}\!\right)}\\ {\dfrac{{p_{i,j,k}^{n + 1} - p_{i,j,k}^{n + \frac{2}{3}}}}{{\Delta t}} = \dfrac{1}{{{h^2}}}\delta _z^2\left( {\dfrac{{p_{i,j,k}^{n + 1} - p_{i,j,k}^n}}{2}} \right)} \end{array}} \right. $

(16)

为了求解压力分布，将油水相基本渗流微分方程(1)和(2)相加并写为三维的形式，推导得到只含有一个未知量${p_{\rm{o}}}$的压力方程

$\begin{align} & {{B}_{\text{o}}}\cdot \text{ }\!\!\{\!\!\text{ }\frac{\partial }{\partial x}\left[ {{\lambda }_{\text{o}x}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial x}-{{F}_{\text{g}}}\cdot {{\rho }_{\text{o}}}\text{g}\frac{\partial D}{\partial x} \right) \right]+\frac{\partial }{\partial y}\left[ \lambda {{}_{\text{o}y}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial y}-{{F}_{\text{g}}}\cdot \rho {{}_{\text{o}}}\text{g}\frac{\partial D}{\partial y} \right) \right] \\ & +\left. \frac{\partial }{\partial z}\left[ \lambda {{}_{\text{o}z}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial z}-{{F}_{\text{g}}}\cdot {{\rho }_{\text{o}}}\text{g}\frac{\partial D}{\partial z} \right) \right]+{{q}_{\text{ov}}} \right\} \\ & {{B}_{\text{w}}}\text{ }\!\!\{\!\!\text{ }\frac{\partial }{\partial x}\left[ \lambda {{}_{\text{w}x}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial x}-{{F}_{\text{g}}}\cdot \rho {{}_{\text{w}}}\text{g}\frac{\partial D}{\partial x} \right) \right]+\frac{\partial }{\partial y}\left[ {{\lambda }_{\text{w}y}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial y}-{{F}_{\text{g}}}\cdot \rho {{}_{\text{w}}}\text{g}\frac{\partial D}{\partial y} \right) \right] \\ & +\frac{\partial }{\partial z}\left[ \lambda {{}_{\text{w}z}}\cdot \left( \frac{\partial {{p}_{\text{o}}}}{\partial z}-{{F}_{\text{g}}}\cdot \rho {{}_{\text{w}}}\text{g}\frac{\partial D}{\partial z} \right) \right]+{{q}_{\text{wv}}}\text{ }\!\!\}\!\!\text{ }=\phi {{C}_{\text{tr}}}\cdot \frac{\partial {{p}_{\text{o}}}}{\partial t} \\ \end{align}$

(17)

流动系数的定义如下

$ \lambda _{{\rm{pd}}} = F_{{\rm{Dc}}} \cdot \dfrac{{K_d K_{{\rm{rp}}} }}{{B_{\rm{p}} \mu _{\rm{p}} }}~~~~~~({\rm{p}} = {\rm{o}}，{\rm{w}}；{\rm{ }}d = x，y，z) $

对压力方程(17)，按照Douglas差分格式式(16)展开即可得到相应的压力线性方程组。

首先定义如下变量

$\begin{align} & {{T}_{\text{o}xi+\frac{1}{2}}}=\frac{\Delta {{y}_{j}}\Delta {{z}_{k}}{{K}_{xi+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{x}_{i+1}}+\Delta {{x}_{i}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{i+\frac{1}{2}}}, \\ & {{T}_{\text{o}xi-\frac{1}{2}}}=\frac{\Delta {{y}_{j}}\Delta {{z}_{k}}{{K}_{xi-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{x}_{i-1}}+\Delta {{x}_{i}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{i-\frac{1}{2}}}, \\ & {{T}_{\text{o}yj+\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{z}_{k}}{{K}_{yj+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{y}_{j+1}}+\Delta {{y}_{j}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{j+\frac{1}{2}}} \\ & {{T}_{\text{o}yj-\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{z}_{k}}{{K}_{yj-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{y}_{j-1}}+\Delta {{y}_{j}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{j-\frac{1}{2}}}, \\ & {{T}_{\text{o}zk+\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{y}_{j}}{{K}_{zk+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{z}_{k+1}}+\Delta {{z}_{k}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{k+\frac{1}{2}}}, \\ & {{T}_{\text{o}zk-\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{y}_{j}}{{K}_{zk-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{z}_{k-1}}+\Delta {{z}_{k}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{ro}}}}{{{B}_{\text{o}}}{{\mu }_{\text{o}}}} \right)}_{k-\frac{1}{2}}} \\ & {{T}_{\text{w}xi+\frac{1}{2}}}=\frac{\Delta {{y}_{j}}\Delta {{z}_{k}}{{K}_{xi+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{x}_{i+1}}+\Delta {{x}_{i}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{i+\frac{1}{2}}}, \\ & {{T}_{\text{w}xi-\frac{1}{2}}}=\frac{\Delta {{y}_{j}}\Delta {{z}_{k}}{{K}_{xi-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{x}_{i-1}}+\Delta {{x}_{i}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{i-\frac{1}{2}}}, \\ & {{T}_{\text{w}yj+\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{z}_{k}}{{K}_{yj+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{y}_{j+1}}+\Delta {{y}_{j}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{j+\frac{1}{2}}} \\ & {{T}_{\text{w}yj-\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{z}_{k}}{{K}_{yj-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{y}_{j-1}}+\Delta {{y}_{j}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{j-\frac{1}{2}}}, \\ & {{T}_{\text{w}zk+\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{y}_{j}}{{K}_{zk+\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{z}_{k+1}}+\Delta {{z}_{k}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{k+\frac{1}{2}}}, \\ & {{T}_{\text{w}zk-\frac{1}{2}}}=\frac{\Delta {{x}_{i}}\Delta {{y}_{j}}{{K}_{zk-\frac{1}{2}}}}{\frac{1}{2}\left( \Delta {{z}_{k-1}}+\Delta {{z}_{k}} \right)}{{\left( \frac{{{F}_{\text{Dc}}}{{K}_{\text{rw}}}}{{{B}_{\text{w}}}{{\mu }_{\text{w}}}} \right)}_{k-\frac{1}{2}}} \\ \end{align}$

对应于Douglas差分格式第一步的三对角形式的压力线性方程组为

$\begin{align} & \frac{1}{2}\left( {{B}_{\text{o}}}{{T}_{\text{o}xi-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi-\frac{1}{2}}} \right)p_{\text{o}i-1,j,k}^{n+\frac{1}{3}}-\frac{1}{2} \\ & \left[ {{B}_{\text{o}}}{{T}_{\text{o}xi+\frac{1}{2}}}+{{B}_{\text{o}}}{{T}_{\text{o}xi-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi-\frac{1}{2}}}+2\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t} \right] \\ & p_{\text{o}i,j,k}^{n+\frac{1}{3}}+\frac{1}{2}\left( {{B}_{\text{o}}}{{T}_{\text{o}xi+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi+\frac{1}{2}}} \right)p_{\text{o}i+1,j,k}^{n+\frac{1}{3}}= \\ & -\frac{1}{2}\left[ \left( {{B}_{\text{o}}}{{T}_{\text{o}xi+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi+\frac{1}{2}}} \right)\left( p_{\text{o}i+1,j,k}^{n}-p_{\text{o}i,j,k}^{n} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{T}_{\text{o}xi-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}xi-\frac{1}{2}}} \right)\left( p_{\text{o}i-1,j,k}^{n}-p_{\text{o}i,j,k}^{n} \right) \right]+ \\ & {{F}_{\text{g}}}\cdot \text{g}\left[ \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}xi+\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}xi+\frac{1}{2}}} \right)\left( {{D}_{i+1,j,k}}-{{D}_{i,j,k}} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}xi-\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}xi-\frac{1}{2}}} \right)\left( {{D}_{i-1,j,k}}-{{D}_{i,j,k}} \right) \right] \\ & -\left[ \left( {{B}_{\text{o}}}{{T}_{\text{o}yj+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj+\frac{1}{2}}} \right)\left( p_{\text{o}i,j+1,k}^{n}-p_{\text{o}i,j,k}^{n} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{T}_{\text{o}yj-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj-\frac{1}{2}}} \right)\left( p_{\text{o}i,j-1,k}^{n}-p_{\text{o}i,j,k}^{n} \right) \right]+ \\ & {{F}_{\text{g}}}\cdot \text{g}\left[ \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}yj+\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}yj+\frac{1}{2}}} \right)\left( {{D}_{i,j+1,k}}-{{D}_{i,j,k}} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}yj-\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}yj-\frac{1}{2}}} \right)\left( {{D}_{i,j-1,k}}-{{D}_{i,j,k}} \right) \right] \\ & -\left[ \left( {{B}_{\text{o}}}{{T}_{\text{o}zk+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk+\frac{1}{2}}} \right)\left( p_{\text{o}i,j,k+1}^{n}-p_{\text{o}i,j,k}^{n} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{T}_{\text{o}zk-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk-\frac{1}{2}}} \right)\left( p_{\text{o}i,j,k-1}^{n}-p_{\text{o}i,j,k}^{n} \right) \right]+ \\ & {{F}_{\text{g}}}\cdot \text{g}\left[ \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}zk+\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}zk+\frac{1}{2}}} \right)\left( {{D}_{i,j,k+1}}-{{D}_{i,j,k}} \right)+ \right. \\ & \left. \left( {{B}_{\text{o}}}{{\rho }_{\text{o}}}{{T}_{\text{o}zk-\frac{1}{2}}}+{{B}_{\text{w}}}{{\rho }_{\text{w}}}{{T}_{\text{w}zk-\frac{1}{2}}} \right)\left( {{D}_{i,j,k-1}}-{{D}_{i,j,k}} \right) \right] \\ & -{{B}_{\text{o}}}{{Q}_{\text{o}i,j,k}}-{{B}_{\text{w}}}{{Q}_{\text{w}i,j,k}}-\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t}p_{\text{o}i,j,k}^{n} \\ \end{align}$

(18)

对应于Douglas差分格式第二步的三对角形式的压力线性方程组为

$\begin{align} & \left( {{B}_{\text{o}}}{{T}_{\text{o}yj-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj-\frac{1}{2}}} \right)p_{\text{o}i,j-1,k}^{n+\frac{2}{3}}- \\ & \left[ {{B}_{\text{o}}}{{T}_{\text{o}yj+\frac{1}{2}}}+{{B}_{\text{o}}}{{T}_{\text{o}yj-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj-\frac{1}{2}}}+2\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t} \right]p_{\text{o}i,j,k}^{n+\frac{2}{3}}+ \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}yj+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj+\frac{1}{2}}} \right)p_{\text{o}i,j+1,k}^{n+\frac{2}{3}}= \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}yj+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj+\frac{1}{2}}} \right)\left( p_{\text{o}i,j+1,k}^{n}-p_{\text{o}i,j,k}^{n} \right)+ \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}yj-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}yj-\frac{1}{2}}} \right)\left( p_{\text{o}i,j-1,k}^{n}-p_{\text{o}i,j,k}^{n} \right)-2\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t}p_{\text{o}i,j,k}^{n+\frac{1}{3}} \\ \end{align}$

(19)

对应于Douglas差分格式第三步的三对角形式的压力线性方程组为

$\begin{align} & \left( {{B}_{\text{o}}}{{T}_{\text{o}zk-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk-\frac{1}{2}}} \right)p_{\text{o}i,j,k-1}^{n+1}- \\ & \left[ {{B}_{\text{o}}}{{T}_{\text{o}zk+\frac{1}{2}}}+{{B}_{\text{o}}}{{T}_{\text{o}yzk-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk-\frac{1}{2}}}+2\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t} \right] \\ & p_{\text{o}i,j,k}^{n+1}+ \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}zk+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk+\frac{1}{2}}} \right)p_{\text{o}i,j,k+1}^{n+1}= \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}zk+\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk+\frac{1}{2}}} \right)\left( p_{\text{o}i,j,k+1}^{n}-p_{\text{o}i,j,k}^{n} \right)+ \\ & \left( {{B}_{\text{o}}}{{T}_{\text{o}zk-\frac{1}{2}}}+{{B}_{\text{w}}}{{T}_{\text{w}zk-\frac{1}{2}}} \right)\left( p_{\text{o}i,j,k-1}^{n}-p_{\text{o}i,j,k}^{n} \right)\quad -2\frac{\left( {{V}_{\text{pr}}}{{C}_{\text{tr}}} \right)_{i,j,k}^{n}}{\Delta t}p_{\text{o}i,j,k}^{n+\frac{2}{3}} \\ \end{align}$

(20)

求解方程组式(18)~式(20)即可得到油相压力分布。在此基础上，将${p_{\rm{o}}}$ 代入毛管压力方程得到水相压力${p_{\rm{w}}}$。对于油、水相饱和度及各组分浓度，可通过油水相及各组分基本微分方程的差分展开式，显式求解得到。

4 实例计算分析

为了验证本文提出的采用Douglas格式差分渗流方程进行模型求解的优越性，选取一个五点法井网的非均质模型进行对比模拟研究，模型大小为385 m×385 m×20 m，划分为11×11×5的均匀网格系统，模型第5层渗透率场分布如图 1所示。模拟时长7 000 d，其中，在2 160~2 700 d注入2 000 mg/L的聚合物和2 000 mg/L的交联剂，注入段塞大小0.15 PV，并保持注采平衡。

本文基于隐式求解压力、显式求解饱和度与组分浓度的方法，采用Douglas差分格式离散化渗流数学模型，通过追赶法求解代数方程组，模拟时间为172 s。同时，在相同生产控制条件下，利用商业化的油藏数值模拟软件CMG对上述模型进行数值模拟运算，CMG软件中，采用自适应隐式方法，代数方程组采用迭代的不完全LU分解法求解，该算例模拟时间为529 s。

图 2给出了利用本文方法和CMG软件计算得到的交联聚合物驱生产动态，可以看出，两者基本吻合。计算实例表明，对于相同的油藏数值模拟问题，不同求解方法的计算效率差别很大。采用Douglas差分格式并结合追赶法能够大幅度提高模型的求解效率。

5 结论

（1）充分考虑交联聚合物驱油机理，建立了考虑因素较为完整的交联聚合物驱数学模型。

（2）基于隐式求解压力、显式求解饱和度与组分浓度的方法，在常规中心差分格式的基础上提出了求解渗流微分方程的Douglas差分格式，利用分步降维的思路求解交联聚合物驱数学模型。

(3) 计算实例表明，相比于商业化的油藏数值模拟软件 CMG，本文采用Douglas差分格式并结合追赶法求解交联聚合物驱数学模型，能够大幅度提高模型的求解效率。

符 号 说 明

${F_{{\rm{Dc}}}}$—— 达西单位换算系数，${F_{{\rm{Dc}}}}$=0.0864；

${K_{\rm{o}}}$，${K_{\rm{w}}}$—— 油、水相渗透率，mD；

${\mu _{\rm{o}}}$，${\mu _{\rm{w}}}$—— 油、水相黏度，mPa$\cdot$s；

t—— 时间步长，s；

${p_{\rm{o}}}$，${p_{\rm{w}}}$—— 油、水相压力，MPa；

${F_{{\rm{g}}}}$—— 重力项单位换算系数，${F_{{\rm{g}}}}$=1.0×10$^{-6}$；

${\rho _{\rm{o}}}$，${\rho _{\rm{w}}}$—— 油、水相密度，kg/m³；

g—— 重力加速度，g=9.8 m/s²；

D— 标高，m，由某一基准面算起的垂直方向深度，与重力加速度方向相同；

${q_{\rm{ov}}}$，${q_{\rm{wv}}}$—— 单位体积岩石在单位时间内注入(或采出)油、水相的地面体积量，1/d；

$\phi$—— 地层岩石孔隙度，无因次；

${S_{\rm{o}}}$，${S_{\rm{w}}}$—— 油藏含油、含水饱和度，无因次；

${B_{\rm{o}}}$，${B_{\rm{w}}}$—— 油、水相体积系数，无因次；

${R_{\rm{Kw}}}$—— 水相渗透率降低系数，无因次；

${C_{l\rm{w}}}$—— 水相中l 组分的浓度，g/m³；

l—— l=P，G，Cr，聚合物、交联聚合物或交联剂组分；

${\beta_{l}}$—— 组分l的扩散系数，m²/d；

${C_{l\rm{ad}}}$—— 组分l的吸附浓度，g/m³；

${R_{l}}$—— 单位时间、单位孔隙介质体积中l组分的反应消耗或生成的质量，g/(m³$\cdot$d)；

${E_{l}}$—— l组分的可及孔隙体积倍数，无因次； ${p_{\rm{cow}}}$—— 油水两相之间的毛管压力，MPa；

${p_{\rm{c}}}$—— 毛管压力，是含水饱和度和相间界面张力的函数，MPa；

${\sigma _{{\rm{wo}}}}$—— 油水两相之间的界面张力，N/m；

${x_{\rm{w}}}$，${x_{\rm{P}}}$，${x_{\rm{G}}}$—— 水相中水、聚合物和交联聚合物组分的摩尔分数，无因次；

$f({x_{\rm{P}}})$，$f({x_{\rm{G}}})$—— 黏度非线性混合函数中聚合物、交联聚合组分对应于摩尔分数${x_{\rm{P}}}$，${x_{\rm{G}}}$的非线性系数，无因次；

${\mu _{\rm{P}}}$，${\mu _{\rm{G}}}$—— 水相中聚合物、交联聚合物组分对应于摩尔分数 ${x_{\rm{P}}}$，${x_{\rm{G}}}$时的黏度，mPa$\cdot$s；

$\mu _{\rm{w}}^{'}$—— 纯水的黏度，mPa$\cdot$s；

a，b—— Langmuir吸附参数，与组分类型、地层水矿化度和地层岩石性质等因素有关，大小可通过实验确定；

${R_{l{\rm{Kw}}}}$—— l组分吸附产生的水相渗透率降低系数，无因次；

${F_{l{\rm{rr}}}}$—— l组分的残余阻力系数，无因次；

${C_{l{\rm{ad}}\max }}$—— l组分的最大吸附浓度，g/m³；

${F_{lJ}}$—— l组分在J方向的扩散通量，g/(m$^3\cdot$d)；

${D_{lJ}}$—— l组分在J方向的扩散系数，m/d；

polymer，clpolymer，Cr$^{3+}$—— 聚合物、交联聚合物和交联剂；

n—— Cr$^{3+}$与聚合物反应的化学当量比，由实验确定；

$\alpha$—— 反应动力学系数，g$^{-3}$$\cdot$m⁹$\cdot$d$^{-1}$；

${C_{\rm{P}}}$，${C_{\rm{G}}}$，${C_{\rm{Cr}}}$—— 聚合物、交联聚合物和交联剂组分的吸附浓度，g/m³；

h—— x，y，z方向的空间步长，m；

$\delta _x^2$，$\delta _y^2$，$\delta _z^2$—— 关于x，y，z的二阶中心差分；

$K_x$，$K_y$，$K_z$—— x，y，z方向的渗透率，mD；

$\lambda _{{\rm{o}}x}$，$\lambda _{{\rm{o}}y}$，$\lambda _{{\rm{o}}z}$—— 油相在x，y，z方向的流动系数；

$\lambda _{{\rm{w}}x}$，$\lambda _{{\rm{w}}y}$，$\lambda _{{\rm{w}}z}$—— 水相在x，y，z方向的流动系数；

${V_{\rm{pr}}}$—— 网格块的孔隙体积，m³；

${C_{\rm{tr}}}$—— 岩石综合压缩系数，MPa$^{-1}$；

${Q_{{\rm{o}}i,j,k}}$，${Q_{{\rm{w}}i,j,k}}$—— 网格块(i ，j，k)在单位时间内注入(或采出)油水相的地面体积量，m³/d。