基于全隐式有限差分的三元复合驱渗流模型数值计算方法

Numerical Computation Method for Solving Flow Model of ASP Flooding Based on Full Implicit Finite Difference

摘要: 针对三元复合驱缺乏科学的数学模型描述驱油机理，数值求解数值稳定性不高、精度差的问题，本文提出了一种基于全隐式有限差分的数值计算方法。首先综合三元驱替剂(碱、表面活性剂和聚合物)对油、水物理化学特性的影响，以及油水渗流方程和吸附扩散方程，建立了三元复合驱数学模型。然后采用全隐式有限差分法，引入块中心网格系统将数学模型从时间和空间上进行离散，从而将原始的偏微分方程组转化为差分代数方程组。最后，采用Newton-Raphson法进行迭代求解，得到系统的状态和输出。为了验证提出方法的效果，针对四注九采三元复合驱实例，分别采用本文方法和CMG数值模拟软件求解，仿真结果表明，本文提出的方法具有较高的求解精度。

Abstract: For alkali-surfactant-polymer (ASP) flooding, there is lack of scientific mathematical model to describe the flooding mechanism, and the numerical computation has bad numerical stability and low accuracy, this paper presents a new numerical computation method based on full implicit finite difference method. Firstly, a comprehensive ASP flooding model is built with considering the influence on the physicochemical characteristic and seepage equation of oil and water and the adsorption diffusion equation which is caused by the adding of displacing agents. Then the full implicit finite difference method is applied to discretize the mechanism in time and space by block center grid system. Further, the original partial differential equations are transformed into a series of difference differential algebraic equations. At last, the problem is solved iteratively by Newton-Raphson method to get the system state and output. To verify the proposed method, a four-injection-nine-production wells example is introduced which is solved by proposed method and CMG software. Simulation result shows that the proposed method has good accuracy.

1. 引言

三次采油技术是指通过向油藏的驱替流体中加入化学剂、混相溶剂、热量等介质，改变油藏流体和岩层的物理和化学性质，使得采出的原油增多的技术，能在二次采油的基础上进一步提高采收率。三次采油技术包括化学驱、混相驱、热力采油、微生物采油等，其中，化学驱又可分为碱驱、表面活性剂驱、聚合物驱、复合驱等 [1] [2] 。传统的三次采油一般利用化学驱、混相驱、热力采油等，但是随着油藏开采的进行，油田地质条件发生改变，采油难度逐年增加，而且开采效果也越来越不理想，于是在众多的三次采油新技术中，三元复合驱脱颖而出，成为最具有应用前景的方法之一 [3] 。

国内外已经有很多关于聚合物驱以及聚合物二元驱、复合驱的数学模型研究。其中，李宜强等 [4] 采用相似原理建立了泡沫复合驱和三元复合驱的相似准则，并利用数值模拟和物理模拟方法对准则进行检验和敏感性分析；Thomas等 [5] 建立了表面活性剂驱油水两相的数学模型；袁士义等 [6] 研究了碱水复合驱的数学模型；侯建等 [7] 对聚合物驱模型及数值模拟方法进行了研究。杨承志等 [1] 指出，碱–表面活性剂–聚合物三元复合驱提高采收率是在碱驱、表面活性剂/聚合物驱和聚合物驱基础上发展起来的。化学驱的数学模型主要包括油水渗流方程、化学驱替剂的吸附扩散方程、油藏及各组分的物化代数方程以及各种约束条件，然而还没有对三元复合驱数学模型展开深入的研究。在求解数学模型时，包括各种数值模拟软件(如VIP，CMG等)，普遍采用有限差分法 [8] ，对偏微分方程进行离散化，转化为代数方程的非线性规划问题进行优化求解，由于问题规模庞大，求解时累积误差严重，极大地影响了结果的准确性。

但随着人们对复合驱驱油机理认识的不断加深，对三元复合驱采油过程提出了更多的问题和要求。一方面现有的复合驱数值模拟软件对三元复合驱的物理化学现象还没有完全描述清楚，对于三元复合驱的偏微分方程模型通常都采用有限差分方法求解，数值稳定性差，且该问题涉及多个复杂高阶偏微分方程，难于求解和应用，且累积误差较大。针对以上问题，本文提出了一种基于全隐式有限差分的三元复合驱渗流模型计算方法。综合考虑三元驱替剂的加入对油藏渗流的影响，建立了全面合理的三元复合驱渗流模型，并采用数值稳定性较好的全隐式有限差分法对数学模型进行离散，转化为差分代数方程组，通过Newton-Raphson法进行求解，得到整个空间状态和输出。

2. 三元复合驱数学模型

三元复合驱的数学模型主要包括油、水渗流方程和驱替剂吸附扩散方程，是一组典型的拋物型偏微分方程组。另外还包括一系列物化代数方程：碱耗方程、渗透率方程、驱替剂吸附方程、粘度方程等。对于实际的矿场，对应的物化代数方程往往根据不同的地质条件，由实验数据确定。

三元复合驱问题的基本假设条件如下：

1) 地层均质，油藏等温，且吸附现象满足Langmuir等温吸附公式；

2) 驱油体系由碱、聚合物、表面活性剂组成，驱油过程为油水两相流，水相中存在碱、聚合物和表面活性剂，油相中存在表面活性剂；

3) Darcy定律适合化学剂存在情况下的油水两相；

4) 各种吸附都达到平衡，相平衡瞬间建立，且满足广义Fick定律；

流体和岩石微可压缩，考虑重力和毛管力的影响，忽略驱替剂对水溶液质量守恒的影响，考虑驱替剂引起的水相渗透率以及水相黏度的变化，并且考虑聚合物的不可及孔隙体积。

基于以上假设条件，设油藏所在的三维区域为 $\Omega \in R^3$ ， $\left(x,y,z\in \Omega \right)$ ， $\partial \Omega$ 为区域的边界， $n$ 表示边界的法线方向，油藏中主要含有n种组分(油、水、聚合物、表面活性剂、碱)，则三元复合驱的数学模型 [1] [3] 如下所示：

油相渗流连续性方程：

$\nabla \cdot \left[\frac{K{k}_{ro}}{B_o{\mu }_o}\nabla \left(p_o-{\rho }_{o}gh\right)\right]+q_o=\frac{\partial }{\partial t}\left[\frac{\varphi \left(1-S_w\right)}{B_o}\right],$ (1)

水相渗流连续性方程：

$\nabla \cdot \left(\frac{K{k}_{rw}}{B_w{R}_k{\mu }_w}\nabla \left(p_w-{\rho }_{w}gh\right)\right)+q_w=\frac{\partial }{\partial t}\left(\frac{\varphi S_w}{B_w}\right),$ (2)

聚合物吸附扩散方程：

$\nabla \cdot \left(\frac{K{k}_{rw}{c}_p}{B_w{R}_k{\mu }_w}\nabla \left(p_w-{\rho }_{w}gh\right)\right)+\nabla \cdot \left[\left(D_w+D_{wp}\right)\frac{{\varphi }_p{S}_w}{B_w}\nabla c_p\right]+q_c=\frac{\partial }{\partial t}\left({\rho }_r\left(1-\varphi \right)C_{rp}+\frac{{\varphi }_p{S}_w{c}_p}{B_w}\right),$ (3)

表面活性剂吸附扩散方程：

$\begin{array}{l}\nabla \cdot \left(\frac{K{k}_{rw}{c}_{ws}}{B_w{R}_k{\mu }_w}\nabla \left(p_w-{\rho }_{w}gh\right)\right)+\nabla \cdot \left(\frac{K{k}_{ro}{c}_{os}}{B_o{R}_k{\mu }_o}\nabla \left(p_o-{\rho }_{o}gh\right)\right)\\ +\nabla \cdot \left[\left(D_o+D_{os}\right)\frac{{\varphi }_s{S}_o}{B_o}\nabla c_{os}\right]+\nabla \cdot \left[\left(D_w+D_{ws}\right)\frac{{\varphi }_s{S}_w}{B_w}\nabla c_{ws}\right]+q_d\\ =\frac{\partial }{\partial t}\left(\frac{{\varphi }_s{S}_o{c}_{os}}{B_o}+\frac{{\varphi }_s{S}_w{c}_{ws}}{B_w}+{\rho }_r\left(1-\varphi \right)C_{rs}\right),\end{array}$ (4)

碱组分吸附扩散方程：

$\begin{array}{l}\nabla \cdot \left(D_{OH}\frac{\varphi S_w}{B_w}\nabla c_{OH}\right)+\nabla \cdot \left(\frac{K{k}_{rw}{c}_{OH}}{B_w{R}_k{\mu }_w}\nabla \left(p_w-{\rho }_{w}gh\right)\right)+R_{OH}+q_e\\ =\frac{\partial }{\partial t}\left[\frac{K_a}{{\left(1+K_b{c}_{OH}\right)}^2}\left(1-\varphi \right)S_w{c}_{OH}+\frac{\varphi S_w{c}_{OH}}{B_w}\right],\end{array}$ (5)

初始条件：

$$ $\begin{array}{l}{p\left(x,y,z,t\right)|}_{t=0}=p^0,{S_w\left(x,y,z,t\right)|}_{t=0}=S_w^0,\\ {c_{\Theta }\left(x,y,z,t\right)|}_{t=0}=c_{\Theta }^0\text{, }\left(x,y,z\right)\in \Omega ,\end{array}$ (6)

边界条件：

${\frac{\partial p}{\partial n}|}_{\partial \Omega }=0,{\frac{\partial S_w}{\partial n}|}_{\partial \Omega }=0,{\frac{\partial c_{\Theta }}{\partial n}|}_{\partial \Omega }=0,$ (7)

上述模型中各符号的物理意义如下：

在所有下标中 $p,s,OH$ 分别表示聚合物、表面活性剂、碱；K为绝对渗透率[μm²]； $k_{ro},k_{rw}$ 分别为油相和水相的相对渗透率[-]； $p_o,p_w$ 为油相、水相压力[MPa]； $S_o,S_w$ 分别为含油饱和度、含水饱和度[-]； $p_{cow}\left(x,y,z,t\right)$ 为毛管力[MPa]；g为重力加速度[m/s²]； $h\left(x,y,z\right)$ 为深度[m]，向下为正； $c_{\Theta }$ 为相应驱替剂的浓度 $\Theta =\left\{p,s,OH\right\}$ 分别代表聚合物浓度[g/L]、表面活性剂质量浓度[%]、碱的质量浓度[%]， $c_{ij}$ 表示j在i中的质量浓度[%]； ${\rho }_o,{\rho }_w,{\rho }_r$ 分别油相、水相和岩石的密度(kg/m³)； $B_o,B_w$ 分别油和水的体积系数[-]； ${\mu }_o,{\mu }_w$ 为加入驱替剂后油相和水相的黏度[mPa×s]； $\varphi ,{\varphi }_p,{\varphi }_s$ 分别为岩石的孔隙度、聚合物可及孔隙度以及表面活性剂可及孔隙度[-]， ${\varphi }_{p,s}=\varphi f_{a,s}$ ， $f_a,f_s$ 为可及孔隙系数[-]； $R_k$ 为相对渗透率下降系数[-]； $K_a$ 为表征离子交换与吸附量大小的参数[-]； $K_b$ 为碱的吸附常数，[cm³/g]； $v_w$ 为渗流速度[cm/s]，流入为正； $R_{OH}$ 为碱耗，[1/d]； $C_{rp},C_{rs}$ 为单位质量岩石吸附聚合物、表面活性剂的质量[mg/g]； $r_o,r_w$ 分别为油相和水相的流动系数[μm²/(mPa×s)]； $q_o$ 为油相在标准状态下的流速[1/d]，流入为正； $q_w$ 为水相在标准状态下的流速[1/d]，流入为正； $q_c,q_d,q_e$ 为井筒驱替剂的运移速度，单位为[g/(d×L)], [1/(100d)], [1/(100d)]； $D_i,i\in \left\{w,o,OH\right\}$ 为扩散系数[m²/s]， $D_{ij},i\in \left\{w,o\right\},j\in \left\{s,p\right\}$ 表示j在i中的扩散系数[m²/s]；t为时间[d]； $\nabla$

为Hamilton算子，在直角坐标系中 $\nabla =\frac{\partial }{\partial x}i+\frac{\partial }{\partial y}j+\frac{\partial }{\partial z}k$ ； $x,y,z$ 为直角坐标系的三个方向，长度单位[m]。

上述代数关系式，可以通过模型中的辅助变量代入支配方程，不需要作为最优控制问题的等式约束。在上述模型中系统的状态变量为 $p,S_w,c_{\Theta }$ ，分别为油藏压力，含水饱和度和油藏中驱替剂浓度(聚合物、表面活性剂和碱)。控制变量(驱替剂的注入浓度)为 $c_{\Theta in}$ ，表示水井中相应驱替剂的注入浓度。 $q_o,q_w$ 分别为油相和水相的流速， $q_c,q_d,q_e$ 为井筒驱替剂的运移速度。 $q_{in}\left(x,y,z,t\right)\ge 0$ 表示在位置 $\left(x,y,z\right)$ 有井的注入项， $q_{out}\left(x,y,z,t\right)\ge 0$ 表示在位置 $\left(x,y,z\right)$ 有井的产出项。 $q_o,q_w,q_c,q_{in},q_{out}$ 以及控制变量 $c_{\Theta in}$ 只在有井的位置有定义。将所有注入项和所有产出项位置的集合分别记为 ${\psi }_w=\left\{{\left(x,y,z\right)}_{j,k},j\in {\kappa }_w,k\in {\vartheta }_w{}_j\right\}$ 和 ${\psi }_p=\left\{{\left(x,y,z\right)}_{j,k},j\in {\kappa }_p,k\in {\vartheta }_p\right\}$ ，则方程中的流量项定义为：

$q_o=\{\begin{array}{l}-\left(1-f_w\right)q_{out},\left(x,y,z\right)\in {\psi }_p\\ 0,\left(x,y,z\right)\notin {\psi }_p\end{array}$ (8)

$q_w=\{\begin{array}{l}-f_w{q}_{out},\left(x,y,z\right)\in {\psi }_p\\ q_{in},\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (9)

$q_c=\{\begin{array}{l}{q}_w{c}_p,\left(x,y,z\right)\in {\psi }_p\\ q_w{c}_{pin},\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (10)

$q_d=\{\begin{array}{l}{q}_w{c}_{ws}+q_o{c}_{os},\left(x,y,z\right)\in {\psi }_p\\ q_w{c}_{sin},\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (11)

$q_e=\{\begin{array}{l}{q}_w{c}_{OH},\left(x,y,z\right)\in {\psi }_p\\ q_w{c}_{OHin},\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (12)

模型中，油、水体积系数均为油藏压力 $p\left(x,y,z,t\right)$ 的函数，具体如下：

$B_o=B_{or}/\left[1+C_o\left(p-p_r\right)\right],$ (13)

$B_w=B_{wr}/\left[1+C_w\left(p-p_r\right)\right],$ (14)

其中， $p_r$ 表示参考压力(MPa)， $B_{wr},B_{or}$ 分别表示参考压力下的水、油体积系数， $C_o,C_w$ 表示油、水压缩系数。

流动系数定义：

$r_o=\frac{K{k}_{ro}}{B_o{\mu }_o},r_w=\frac{K{k}_{rw}}{B_w{R}_k{\mu }_w}$ (15)

$f_w\left(x,y,z,t\right)$ 为含水率是状态 $p,S_w,c_p,c_s,c_{OH}$ 的函数，其定义为 [3]

$f_w=\frac{r_w}{r_w+r_o},$ (16)

表面活性剂浓度：

$c_s=\frac{q_w{c}_{ws}+q_o{c}_{os}}{q_w+q_o},$ (17)

物化代数方程详见附录。

3. 全隐式有限差分求解

三元复合驱的机理模型(1)~(5)是一组复杂的非线性多变量耦合偏微分方程组，其中，状态变量分别为油藏压力 $p\left(x,y,z,t\right)$ 、含水饱和度 $S_w\left(x,y,z,t\right)$ 以及驱替剂溶液浓度 $c_{\Theta }\left(x,y,z,t\right),\Theta =\left\{OH,s,p\right\}$ 。目前主要采用有限差分法求解偏微分方程组，该方法主要包括两类：一类是依次求解各个状态变量，如隐压显饱法，该方法先隐式求解状态压力，再显式求解含水饱和度、驱替剂浓度；另一类是联立求解方程组，同时求解各个状态变量，如全隐式法。隐压显饱法虽然在每个时间步需要的计算量较小，但是为了保证求解的稳定性，步长必须取得很小，有时候甚至小到没有实际价值。全隐式法离散化模型后得到的差分方程组往往阶数较高，需要的计算量较大，但是该方法可以选取的步长较大，具有很好的数值稳定性。此外，该方法离散三元复合驱数学模型得到的最优控制模型和伴随方程具有特殊的形式，便于编程求解。随着计算机技术的飞速发展，计算机的运行速度和存储能力都有了极大提高，结合矩阵存储以及数据处理等数值求解技术，完全可以实现全隐式方法的计算。因此，本文采用全隐式有限差分法 [9] 求解三元复合驱数学模型。

3.1. 全隐式有限差分离散化

对于(1)~(5)所述的三元复合驱数学模型，优化变量为注入井的驱替剂(碱、表面活性剂和聚合物)注入浓度 $c_{\Theta in}$ ，状态变量为油藏空间每一个空间点的压力 $p\left(x,y,z,t\right)$ 、含水饱和度 $S_w\left(x,y,z,t\right)$ 和驱替剂的网格浓度 $c_{\Theta }\left(x,y,z,t\right)$ ，系统的输出为采出井的含水率。本节中，我们采用全隐式有限差分法进行离散 [9] ，为了使本问题的数学描述同一般优化问题一致，引入 $x={\left[c_{\Theta },s_w,p\right]}^{\text{T}}$ 表示所有的空间状态变量， $u=c_{\Theta in},\Theta =\left\{p,s,OH\right\}$ 表示所有的优化变量。

将整个油藏离散化为 $n_x\times n_y\times n_z$ 个网格，其中整数 $i\left(i=1,2,\cdots ,n_x\right)$ ， $j\left(j=1,2,\cdots ,n_y\right)$ ， $k\left(k=1,2,\cdots ,n_z\right)$ 分别表示在 $x,y,z$ 方向上的网格序号，即 $x_i,y_j,z_k$ 分别表示对应的 $x,y,z$ 方向上的第 $i,j,k$ 个坐标值。按照块中心网格系统 [10] 进行离散化，即网格点 $\left(x_i,y_j,z_k\right)$ 位于网格 $\left(i,j,k\right)$ 的中心。那么，网格点 $\left(x_i,y_j,z_k\right)$ 处的离散状态向量可表述为

$x_{i,j,k}=x\left(x_i,y_j,z_k\right)\text{.}$ (18)

对于网格 $\left(x_i,y_j,z_k\right)$ ，定义 $x_{i-1/2},x_{i+1/2}$ ， $y_{j-1/2},y_{j+1/2}$ ， $z_{k-1/2},z_{k+1/2}$ 分别为 $x,y,z$ 方向上的左、右边界坐标。仅在状态网格上离散化三元复合驱机理模型，根据网格系统的划分结果，状态变量的有限差分可以表述为

$\frac{\partial x}{\partial t}=\frac{x^{n+1}-x^n}{\Delta t^n},$ (19)

$\frac{\partial x}{\partial x}=\frac{x_{i+1/2,j,k}-x_{i-1/2,j,k}}{\Delta x},$ (20)

$\frac{\partial x}{\partial y}=\frac{x_{i,j+1/2,k}-x_{i,j-1/2,k}}{\Delta y},$ (21)

$\frac{\partial x}{\partial z}=\frac{x_{i,j,k+1/2}-x_{i,j,k-1/2}}{\Delta z},$ (22)

其中， $n=0,1,\cdots ,N-1$ 表示三元复合物的时间步，N表示总的步数， $\Delta t^n$ 表示步数n的长度， $x^{n+1}$ 为时间 $\Delta t^{n+1}$ 处的状态向量， $\Delta t^{n+1}$ 表示第n个时间步终点时刻的仿真时间( $t^N=t_f$ )。 $\Delta x,\Delta y,\Delta z$ 分别表示在 $x,y,z$ 方向上的空间步长。网格在 $x-y$ 平面上的详细划分情况如图1所示。

Figure 1. Schematic diagram of grid division for finite difference

图1. 有限差分网格划分原理图

采用上述有限差分法对三元复合驱数学模型(1)~(5)进行离散化，并在方程组的每一个等式两端同时乘以网格体积 $V=\Delta x\Delta y\Delta z\left({\text{m}}^3\right)$ 。采用全隐式有限差分法将模型中二阶导数项离散化后，状态变量均取 $\Delta t^{n+1}$ 时刻的值。则初始的三元复合驱数学模型可以被离散化为如下一般形式

$G^n=F^{n+1}\left({\tilde{x}}^{n+1}\right)+W^{n+1}\left({\tilde{x}}^{n+1},u^n\right)-\left[A^{n+1}\left({\tilde{x}}^{n+1}\right)-A^n\left({\tilde{x}}^n\right)\right]=0,$ (23)

其中， $G={\left[G_{1,1,1}^{\text{T}},\cdots ,G_{i,j,k}^{\text{T}},\cdots ,G_{n_x,n_y,n_z}^{\text{T}}\right]}^{\text{T}}\in R^{n_x\times n_y\times n_z\times 3}$ 表示离散的数学模型差分方程组，

$\tilde{x}={\left[{\tilde{x}}_{1,1,1}^{\text{T}},\cdots ,{\tilde{x}}_{i,j,k}^{\text{T}},\cdots ,{\tilde{x}}_{n_x,n_y,n_z}^{\text{T}}\right]}^{\text{T}}\in R^{n_x\times n_y\times n_z\times 3}$ 表示离散的状态向量， $u={\left[u_{i,j}|u_{i,j}=c_{\Theta in,i,j},\left(i,j\right)\in {\kappa }_w\right]}^{\text{T}}\in R^{N_w}$ 表示所有注入井的驱替剂注入浓度， ${\kappa }_w$ 是注入井的位置集合， $F$ 表示差分方程组的流动项， $W$ 表示源汇项， $A$ 表示累积项。

在 $t^n$ 时刻网格 $\left(i,j,k\right)$ 的差分方程组 $G_{i,j,k}={\left[G_{oi,j,k},G_{wi,j,k},G_{pi,j,k},G_{si,j,k},G_{OHi,j,k}\right]}^{\text{T}}$ 为

$G_{\Psi }{}_{i,j,k}^n=F_{\Psi }{}_{i,j,k}^{n+1}+W_{\Psi }{}_{i,j,k}^{n+1}-A_{\Psi }{}_{i,j,k}^{n+1}+A_{\Psi }{}_{i,j,k}^n=0,\Psi =\left\{w,o,OH,s,p\right\},$ (24)

该方程实际是五个方程的综合， $\Psi$ 中的 $w,o,OH,s,p$ 分别对应离散水相方程、油相方程、以及碱、表面活性剂和聚合物的吸附扩散方程。

根据式(15)以及第(8)~(12)中的流动项系数 $r_{{\rm Z}},{\rm Z}=\left\{o,w\right\},r_{c\Theta },r_{d\Theta },\Theta =\left\{OH,s,p\right\}$ 定义，流动项 $F_{i,j,k}=[F_{oi,j,k},F_{wi,j,k},F_{pi,j,k},{F_{si,j,k},F_{OHi,j,k}]}^{\text{T}}$ 的全隐式有限差分可写作

$\begin{array}{c}{F}_{{\rm Z}}{}_{i,j,k}^{n+1}=\Delta y\Delta z\left[\frac{r_{{\rm Z}}{}_{i+\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(E_{i+1,j,k}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{{\rm Z}}{}_{i-\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(E_{i-1,j,k}^{n+1}-E_{i,j,k}^{n+1}\right)\right]\\ +\Delta x\Delta z\left[\frac{r_{{\rm Z}}{}_{i,j+\frac{1}{2},k}^{n+1}}{\Delta y}\left(E_{i,j+1,k}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{{\rm Z}}{}_{i,j-\frac{1}{2},k}^{n+1}}{\Delta y}\left(E_{i,j-1,k}^{n+1}-E_{i,j,k}^{n+1}\right)\right]\\ +\Delta x\Delta y\left[\frac{r_{{\rm Z}}{}_{i,j,k+\frac{1}{2}}^{n+1}}{\Delta z}\left(E_{i,j,k+1}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{{\rm Z}}{}_{i,j,k-\frac{1}{2}}^{n+1}}{\Delta z}\left(E_{i,j,k-1}^{n+1}-E_{i,j,k}^{n+1}\right)\right],\end{array}$ (25)

$\begin{array}{c}{F}_{\Theta }{}_{i,j,k}^{n+1}=\Delta y\Delta z\left[\frac{r_{d\Theta }{}_{i+\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(c_{\Theta }{}_{i+1,j,k}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)+\frac{r_{d\Theta }{}_{i-\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(c_{\Theta }{}_{i-1,j,k}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)\right]\\ +\Delta x\Delta z\left[\frac{r_{d\Theta }{}_{i,j+\frac{1}{2},k}^{n+1}}{\Delta y}\left(c_{\Theta }{}_{i,j+1,k}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)+\frac{r_{d\Theta }{}_{i,j-\frac{1}{2},k}^{n+1}}{\Delta y}\left(c_{\Theta }{}_{i,j-1,k}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)\right]\\ +\Delta x\Delta y\left[\frac{r_{d\Theta }{}_{i,j,k+\frac{1}{2}}^{n+1}}{\Delta z}\left(c_{\Theta }{}_{i,j,k+1}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)+\frac{r_{d\Theta }{}_{i,j,k-\frac{1}{2}}^{n+1}}{\Delta z}\left(c_{\Theta }{}_{i,j,k-1}^{n+1}-c_{\Theta }{}_{i,j,k}^{n+1}\right)\right]\\ +\Delta y\Delta z\left[\frac{r_{c\Theta }{}_{i+\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(E_{i+1,j,k}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{c\Theta }{}_{i-\frac{1}{2},j,k}^{n+1}}{\Delta x}\left(E_{i-1,j,k}^{n+1}-E_{i,j,k}^{n+1}\right)\right]\end{array}$

$\begin{array}{c}+\Delta x\Delta z\left[\frac{r_{c\Theta }{}_{i,j+\frac{1}{2},k}^{n+1}}{\Delta y}\left(E_{i,j+1,k}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{c\Theta }{}_{i,j-\frac{1}{2},k}^{n+1}}{\Delta y}\left(E_{i,j-1,k}^{n+1}-E_{i,j,k}^{n+1}\right)\right]\\ +\Delta x\Delta y\left[\frac{r_{c\Theta }{}_{i,j,k+\frac{1}{2}}^{n+1}}{\Delta z}\left(E_{i,j,k+1}^{n+1}-E_{i,j,k}^{n+1}\right)+\frac{r_{c\Theta }{}_{i,j,k-\frac{1}{2}}^{n+1}}{\Delta z}\left(E_{i,j,k-1}^{n+1}-E_{i,j,k}^{n+1}\right)\right].\end{array}$ (26)

其中， $E=p-\rho gh$ 表示流体的流势。

结合公式(8)~(12)以及(17)，源汇项 $W_{i,j,k}={\left[W_{oi,j,k},W_{wi,j,k},W_{pi,j,k},W_{si,j,k},W_{OHi,j,k}\right]}^{\text{T}}$ 可表示为

$W_o{}_{i,j,k}^{n+1}=q_o{}_{i,j,k}^{n+1}=\{\begin{array}{l}-\left(1-f_w{}_{i,j,k}^{n+1}\right)q_{out}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_p\\ 0,\left(x,y,z\right)\notin {\psi }_p\end{array}$ (27)

$W_w{}_{i,j,k}^{n+1}=q_w{}_{i,j,k}^{n+1}=\{\begin{array}{l}-f_w{}_{i,j,k}^{n+1}{q}_{out}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_p\\ q_{in}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (28)

$W_p{}_{i,j,k}^{n+1}=q_c{}_{i,j,k}^{n+1}=\{\begin{array}{l}-f_w{}_{i,j,k}^{n+1}{q}_{out}{}_{i,j}^n{c}_p{}_{i,j,k}^{n+1},\left(x,y,z\right)\in {\psi }_p\\ q_{in}{}_{i,j}^n{c}_{pin}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (29)

$W_s{}_{i,j,k}^{n+1}=q_d{}_{i,j,k}^{n+1}=\{\begin{array}{l}-f_w{}_{i,j,k}^{n+1}{q}_{out}{}_{i,j}^n{c}_{ws}{}_{i,j,k}^{n+1}-\left(1-f_w{}_{i,j,k}^{n+1}\right)q_{out}{}_{i,j}^n{c}_{os}{}_{i,j,k}^{n+1},\left(x,y,z\right)\in {\psi }_p\\ q_{in}{}_{i,j}^n{c}_{sin}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (30)

$W_{OH}{}_{i,j,k}^{n+1}=R_{OH}{}_{i,j,k}^{n+1}+q_e{}_{i,j,k}^{n+1}=\{\begin{array}{l}{R}_{OH}{}_{i,j,k}^{n+1}-f_w{}_{i,j,k}^{n+1}{q}_{out}{}_{i,j}^n{c}_{OH}{}_{i,j,k}^{n+1},\left(x,y,z\right)\in {\psi }_p\\ R_{OH}{}_{i,j,k}^{n+1}+q_{in}{}_{i,j}^n{c}_{OHin}{}_{i,j}^n,\left(x,y,z\right)\in {\psi }_w\\ 0,\left(x,y,z\right)\notin {\psi }_p\cup {\psi }_w\end{array}$ (31)

其中，式(31)为修正碱耗后的碱吸附方程源汇项， $R_{OH}$ 为系统碱耗。结合井参数，系统的含水率 $f_w{}_{i,j,k}^{n+1}$ 可修正为

$f_w{}_{i,j,k}^{n+1}=\frac{W{I}_{i,j,k}\left(r_{wi,j,k}+r_{oi,j,k}\right)}{{\displaystyle \underset{k^{\prime }=1}{\overset{n_z}{\sum }}W{I}_{i,j,k^{\prime }}\left(r_{wi,j,k^{\prime }}+r_{oi,j,k^{\prime }}\right)}},$ (32)

WI表示井指数 [11] ，定义如下

$WI=\frac{2\text{π}K\Delta z}{\ln \frac{d_e}{d_w}+S},$ (33)

其中， $d_e,d_w$ 分别为井筒半径和井筒等效半径，S为表皮系数，当采用矩形离散化网格且渗透率非均质时 $d_e$ 表示如下 [12]

$d_e\cong 0.28{\left[{\left(\frac{K_y}{K_x}\right)}^{\frac{1}{2}}\Delta x^2+{\left(\frac{K_x}{K_y}\right)}^{\frac{1}{2}}\Delta y^2\right]}^{\frac{1}{2}}/{\left(\frac{K_y}{K_x}\right)}^{\frac{1}{2}}+{\left(\frac{K_x}{K_y}\right)}^{\frac{1}{2}},$ (34)

其中， $K_x,K_y$ 分别表示在 $x,y$ 方向上的绝对渗透率。

类似地，采用全隐式有限差分可将累积项 $A_{i,j,k}=[A_{oi,j,k},A_{wi,j,k},A_{pi,j,k},{A_{si,j,k},A_{OHi,j,k}]}^{\text{T}}$ 表示为

$A_o{}_{i,j,k}^{n+1}=\frac{V_{i,j,k}}{\Delta t^n}\frac{{\varphi }_{i,j,k}^{n+1}\left(1-S_w{}_{i,j,k}^{n+1}\right)}{B_o{}_{i,j,k}^{n+1}},A_o{}_{i,j,k}^n=\frac{V_{i,j,k}}{\Delta t^n}\frac{{\varphi }_{i,j,k}^n\left(1-S_w{}_{i,j,k}^n\right)}{B_o{}_{i,j,k}^n},$ (35)

$A_w{}_{i,j,k}^{n+1}=\frac{V_{i,j,k}}{\Delta t^n}\frac{{\varphi }_{i,j,k}^{n+1}{S}_w{}_{i,j,k}^{n+1}}{B_w{}_{i,j,k}^{n+1}},A_w{}_{i,j,k}^n=\frac{V_{i,j,k}}{\Delta t^n}\frac{{\varphi }_{i,j,k}^n{S}_w{}_{i,j,k}^n}{B_w{}_{i,j,k}^n},$ (36)

$\begin{array}{l}{A}_p{}_{i,j,k}^{n+1}=\frac{V_{i,j,k}}{\Delta t^n}\left[{\rho }_r\left(1-{\varphi }_{i,j,k}^{n+1}\right)C_{rp}{}_{i,j,k}^{n+1}+\frac{{\varphi }_p{}_{i,j,k}^{n+1}{S}_w{}_{i,j,k}^{n+1}{c}_p{}_{i,j,k}^{n+1}}{B_w{}_{i,j,k}^{n+1}}\right],\\ A_p{}_{i,j,k}^n=\frac{V_{i,j,k}}{\Delta t^n}\left[{\rho }_r\left(1-{\varphi }_{i,j,k}^n\right)C_{rp}{}_{i,j,k}^n+\frac{{\varphi }_p{}_{i,j,k}^n{S}_w{}_{i,j,k}^n{c}_p{}_{i,j,k}^n}{B_w{}_{i,j,k}^n}\right],\end{array}$ (37)

$\begin{array}{l}{A}_s{}_{i,j,k}^{n+1}=\frac{V_{i,j,k}}{\Delta t^n}\left[\frac{{\varphi }_s{}_{i,j,k}^{n+1}{S}_o{}_{i,j,k}^{n+1}{c}_{os}{}_{i,j,k}^{n+1}}{B_o{}_{i,j,k}^{n+1}}+\frac{{\varphi }_s{}_{i,j,k}^{n+1}{S}_w{}_{i,j,k}^{n+1}{c}_{ws}{}_{i,j,k}^{n+1}}{B_w{}_{i,j,k}^{n+1}}+{\rho }_r\left(1-{\varphi }_{i,j,k}^{n+1}\right)C_{rs}{}_{i,j,k}^{n+1}\right],\\ A_s{}_{i,j,k}^n=\frac{V_{i,j,k}}{\Delta t^n}\left[\frac{{\varphi }_s{}_{i,j,k}^n{S}_o{}_{i,j,k}^n{c}_{os}{}_{i,j,k}^n}{B_o{}_{i,j,k}^n}+\frac{{\varphi }_s{}_{i,j,k}^n{S}_w{}_{i,j,k}^n{c}_{ws}{}_{i,j,k}^n}{B_w{}_{i,j,k}^n}+{\rho }_r\left(1-{\varphi }_{i,j,k}^n\right)C_{rs}{}_{i,j,k}^n\right],\end{array}$ (38)

$\begin{array}{l}{A}_{OH}{}_{i,j,k}^{n+1}=\frac{V_{i,j,k}}{\Delta t^n}\left[\frac{K_a}{{\left(1+K_b{c}_{OH}{}_{i,j,k}^{n+1}\right)}^2}\left(1-{\varphi }_{i,j,k}^{n+1}\right)S_w{}_{i,j,k}^{n+1}{c}_{OH}{}_{i,j,k}^{n+1}+\frac{{\varphi }_{i,j,k}^{n+1}{S}_w{}_{i,j,k}^{n+1}{c}_{OH}{}_{i,j,k}^{n+1}}{B_w{}_{i,j,k}^{n+1}}\right],\\ A_{OH}{}_{i,j,k}^n=\frac{V_{i,j,k}}{\Delta t^n}\left[\frac{K_a}{{\left(1+K_b{c}_{OH}{}_{i,j,k}^n\right)}^2}\left(1-{\varphi }_{i,j,k}^n\right)S_w{}_{i,j,k}^n{c}_{OH}{}_{i,j,k}^n+\frac{{\varphi }_{i,j,k}^n{S}_w{}_{i,j,k}^n{c}_{OH}{}_{i,j,k}^n}{B_w{}_{i,j,k}^n}\right],\end{array}$ (39)

边界条件公式(7)可以表示为

$\frac{\partial x}{\partial x}=0,\frac{\partial x}{\partial y}=0,\frac{\partial x}{\partial z}=0,\forall \left(x,y,z\right)\in \partial \Omega .$ (40)

上式可采用全隐式有限差分离散为

$x_{0,j,k}=x_{1,j,k},x_{n_x+1,j,k}=x_{n_x,j,k},j=1,2,\cdots ,n_y,k=1,2,\cdots ,n_z,$ (41)

$x_{i,0,k}=x_{i,1,k},x_{i,n_y+1,k}=x_{i,n_y,k},i=1,2,\cdots ,n_x,k=1,2,\cdots ,n_z,$ (42)

$x_{i,j,0}=x_{i,j,1},x_{i,j,n_z+1}=x_{i,j,n_z},i=1,2,\cdots ,n_x,j=1,2,\cdots ,n_y.$ (43)

初始条件公式(6)可用全隐式有限差分离散化为

$x_{i,j,k}^0=x\left(x_i,y_j,z_k,0\right),i=1,\cdots ,n_x,j=1,\cdots ,n_y,k=1,\cdots ,n_z.$ (44)

至此，在时域和空间域上连续的三元复合驱偏微分系统公式(1)~(7)已经通过全隐式有限差分法离散化为一系列非线性隐式差分代数方程组。

3.2. 数学模型方程组求解

针对上节得到的差分代数方程组，这里采用Newton-Raphson法 [13] 进行求解：

${\frac{\partial G^n}{\partial {\tilde{x}}^{n+1}}|}_{{\tilde{x}}^{n+1,l}}\delta {\tilde{x}}^{n+1,l}=-G^n\left({\tilde{x}}^{n+1},{\tilde{x}}^n,u^n\right),$ (45)

${\tilde{x}}^{n+1,l+1}={\tilde{x}}^{n+1,l}+\delta {\tilde{x}}^{n+1,l},$ (46)

式中，l表示迭代次数， ${\tilde{x}}^{n+1,l}$ 表示第l次迭代的离散状态向量， ${\frac{\partial G^n}{\partial {\tilde{x}}^{n+1}}|}_{{\tilde{x}}^{n+1,l}}$ 为Jacobian矩阵，计算公式如下：

${\frac{\partial G^n}{\partial {\tilde{x}}^{n+1}}|}_{{\tilde{x}}^{n+1,l}}=\frac{\partial F^{n+1}}{\partial {\tilde{x}}^{n+1,l}}+\frac{\partial W^{n+1}}{\partial {\tilde{x}}^{n+1,l}}-\frac{\partial A^{n+1}}{\partial {\tilde{x}}^{n+1,l}}.$ (47)

通过求解公式(45)可得 $\delta {\tilde{x}}^{n+1,l}$ ， ${\tilde{x}}^{n+1,l+1}$ 的更新通过公式(46)实现。算法每一步迭代的具体实现过程如图2所示，其中，ε为设定的精度。在得到的隐式差分代数方程组的基础上，通过迭代求解，得到整个空间的状态。

4. 仿真实例

4.1. 油藏描述

对于第2节中给出的三维三元复合驱模型，假设三元复合驱油藏由四口注入井和九口采出井构成，所有的井位均匀分布，在每四口采出井中心有一口注入井。其中S表示采出井，I表示注入井。具体井位分布同图3所示。三维油藏参数：长630 m，宽630 m，厚19.990 m，总共有七层，每层厚2.857 m，净厚度为1.4286 m，油层上表面距地表2420 m，每层的孔隙度为0.3，孔隙体积为 $1.1097\times {10}^6{\text{m}}^3$ 。油藏中三种驱替剂(碱、表面活性剂和聚合物)的初始网格浓度为0 g/L。初始渗透率、初始压力和初始含水饱和

Figure 2. Solution process of Newton-Raphson method

图2. Newton-Raphson迭代求解原理图

Figure 3. Distribution of well position

图3. 井位分布图

度如图4、图5和图6所示。整个油藏分为 $x,y,z$ 三个方向，其中 $x,y$ 方向各等分为21个网格，z方向7个网格，总网格数为 $21\times 21\times 7$ 。每口注入井的注入率 $q_{in}$ 为83m³/天，采出井的采出率 $q_{out}$ 如表1所示，部分油藏参数见表2。整个驱油过程持续96个月，三元复合驱周期为48个月，后面为水驱，考虑到三元复合驱注入的初始时间点，总共有 $L=97$ 个时间点，总空间点为 $N=21\times 21\times 7$ 。

4.2. 仿真结果

采用本文提出的方法对三元复合驱模型进行求解，为四口注入井设置相同的注入策略，对于碱、表面活性剂和聚合物的具体注入浓度如下：

$c_{OHin}=\left(3.7,2.8,1.3\right)\text{g}/\text{L},c_{sin}=\left(2.9,1.3,0.6\right)\text{g}/\text{L},c_{pin}=\left(2.7,1.4,1.1\right)\text{g}/\text{L}.$

整个采油周期为96个月，驱替剂注入为三段塞，每个段塞持续16个月，共48个月，剩余时间为水驱。具体的注入策略如图7(a)所示，由于三维三元复合驱中考虑时间、空间和含水饱和度因素，总共五个维度，无法再三维图形中显示。为了说明建模效果，我们选取第三层中从 $\left(1,11,3\right)$ 到 $\left(21,11,3\right)$ 的21个网格，将结果显示在图7中。图7(b)为CMG数值模拟软件得到的空间含水饱和度，图7(c)为采用本文提出方法求解得到的含水饱和度。图8为含水饱和度的求解误差。图9为九口采出井的含水率结果对比。通过对比可以发现，本文提出的方法的求解结果与CMG2010软件的求解结果，无论是对于含水饱和度状态还是含水率都基本一致。从而说明了本文提出方法的有效性。

5. 结论

三元复合驱是一项重要的三次采油技术，由于其模型复杂，涉及变量多，油藏成分不确定等因素，目前很少有学者对三元复合驱的数学模型和计算方法进行系统的研究。本文建立了三元复合驱的渗流模

Figure 4. Distribution of initial permeability

图4. 初始渗透率分布

Figure 5. Distribution of initial pressure

图5. 初始压力分布

Figure 6. Distribution of initial water saturation

图6. 初始含水饱和度分布

(a) (b) (c)

Figure 7. Comparison of solution results. (a) Injection concentration; (b) Simulation result of CMG software; (c) Solution reult with proposed method in this paper

图7. 求解结果对比。(a) 注入浓度；(b) CMG软件结果；(c) 本文提出方法的求解结果

Figure 8. Modeling error of water saturation

图8. 含水饱和度建模误差

Figure 9. Comparison of moisture content for nine production wells

图9. 九口采出井含水率比较

Table 1. Liquid withdrawal amount of production wells

表1. 采出井采液量

Table 2. Partial reservoir parameters for ASP flooding

表2. 三元复合驱部分油藏参数

注：变量的单位为： $M_{\Theta \max }$ -kg， $c_{\Theta \max }$ -g/L，ρ-kg/m³，D-m²/s，μ-μm²，S-[-]。

型，能够科学的描述油水渗流机理，碱、表面活性剂和聚合物的吸附过程，以及对油藏物理化学特性的影响。采用全隐式有限差分法多数学模型进行离散求解，具有较好的数值稳定性和收敛性。将偏微分方程系统离散化为差分方程组后，结合Newton-Raphson法进行求解，通过与成熟数值模拟软件CMG的结果对比，验证了本文提出方法的精度和求解效果。

资助信息

国家自然科学基金(61573378)；山东省自然科学基金(ZR2011FM002)；中央高校基本科研业务费专项基金(15CX06064A)。

附录

与传统的单一化学驱相比，三元复合驱的复杂性不仅表现在支配方程上，多种驱替剂的相互作用，使得物化代数方程也需进行全面的修正。

1) 碱耗

碱耗是三元复合驱中一个非常重要的因素，它直接影响驱油的成败。影响碱耗的因素较多，且较复杂，本模型只考虑驱油中对碱耗影响重要的因素(吸附、酸性物质消耗、重金属离子消耗)，进行简化处理，以达到既能较好的描述碱耗，又能具有较强的实用性的目的。模型中的碱耗通过化学反应项 $r_i$ 进行描述。表达如下 [1] ：

$R_{OH}=-\varphi S_w\frac{\partial }{\partial t}\left(r_1+r_2+r_3\right),$ (48)

式中，r为单位体积的碱消耗量。

考虑的主要影响如下：

a) 溶液中的Na⁺与岩石表面的 ${\text{H}}^{\text{+}}$ 得离子交换引起的快速碱耗 $r_1$ ：该碱耗可近似表示为类似Langmuir型的吸附等温式：

$r_1=r_1{}^0\frac{a_1{c}_{OH}}{1+a_1{c}_{OH}},$ (49)

式中 $r_1{}^0$ 表示该现象引起的最大碱耗， $a_1$ 为系数，均由实验资料确定。

b) 原油中酸性物质引起的碱耗 $r_2$ ：

$r_2=r_2\left({\text{HA}}_w,c_{OH}\right),$ (50)

实验给出不同酸、碱浓度下的碱耗曲线。

c) Ca²⁺，Mg²⁺离子引起的碱耗 $r_3$ ：

$r_3=r_3\left(K_{{\text{Ca}}^{\text{2+}}},K_{{\text{Mg}}^{\text{2+}}},c_{OH}\right),$ (51)

式中 $K_{{\text{Ca}}^{\text{2+}}},K_{{\text{Mg}}^{\text{2+}}}$ 分别表示相应物质的溶度积。

2) 表面张力

采用表面张力等值图描述多种化学剂的复合协同效应，对于给定的原油和配置水的矿化度，表面张力主要随碱、表面活性剂、聚合物的浓度变化 [14] ：

$\sigma =\sigma \left(c_{OH},c_p,c_s\right).$ (52)

3) 毛细管压力

复合体系的毛细管压力可以描述成油水毛细管压力和界面张力的函数：

$p_{cow}\left(S_w\right)=C_{PC}\cdot \sqrt{\frac{\varphi }{K_a}}\cdot {\left(1-S_n\right)}^{N_{PC}},$ (53)

式中 $C_{PC},N_{PC}$ 为常数， $S_n$ 为湿相饱和度。

4) 表面活性剂吸附

表面活性剂的吸附滞留损失可用下式进行描述 [15] ：

$C_{rs}=C_{rs}^0\frac{a_s{c}_s}{1+a_s{c}_s}\left(1-b_s\frac{\text{pH}-\text{7}}{{\text{pH}}_{\text{max}}-\text{7}}\right),$ (54)

式中 $C_{rs}^0,a_s$ 与离子强度有关， ${\text{pH}}_{\text{max}}$ 为注入碱浓度的 $\text{pH}$ 值， $b_s$ 为系数。

5) 聚合物的吸附量

遵循Langmuir等温吸附规律，采用下式进行计算，考虑到聚合物的吸附与盐度的关系为可逆的，与浓度的关系为不可逆的 [16] 。

$C_{rp}=C_{rp}^{\text{max}}\frac{a_1{c}_p}{1+b_1{c}_p},$ (55)

式中， $C_{rp}^{\text{max}}$ 表示不同含盐度下聚合物在岩石表面的最大吸附量， $a_1,b_1$ 为平衡吸附常数。

6) 相对渗透率

水、油相对渗透率 $k_{rw},k_{ro}$ 通常采用插值的方法获得，但是由于插值区间和精度的限制，该方法往往无法获取整个区间的精确值。这里采用辨识的方法获得，具体如下 [17] ：

$K_{rw,ro}=A\cdot {\left(1-S_w\right)}^B\cdot {\left(S_w-C\right)}^D,$ (56)

其中， $A,B,C,D$ 表示通过样本数据辨识得到的参数。

由于多元驱替剂的加入，必须要考虑渗透率的衰减，它主要是由聚合物的吸附滞留所引起，可利用下式进行描述：

$R_k=1+\frac{\left(R_k^{\max }-1\right)\cdot q_p}{q_p^{\max }},$ (57)

式中 $q_p,R_k$ 分别表示不同含盐度下聚合物吸附滞留量和水相渗透率下降系数。

7) 液相的黏度

液相的静止黏度(零剪切速率下的黏度)是聚合物的浓度和含盐度的函数，Flory-Huggins方程用于考虑含盐度的变化情况。

${\mu }_p^0={\mu }_w\left(1+a_1{c}_p^1+a_2{c}_p^2+a_3{c}_p^3+\cdots \right),$ (58)

式中 $a_1,a_2,a_3,\cdots$ 为经验常数，与含盐度有关。

ASP溶液的黏度与剪切速率( $\stackrel{\dot{}}{\gamma }$ )的关系可由Meter方程表达为

${\mu }_{\text{ASP}}={\mu }_w+\frac{{\mu }_p^0-{\mu }_w}{1+{\left[\frac{\stackrel{\dot{}}{\gamma }}{{\stackrel{\dot{}}{\gamma }}_{1/2}}\right]}^{P_a-1}},$ (59)

式中， ${\stackrel{\dot{}}{\gamma }}_{1/2}$ 为当黏度等于 ${\mu }_p^0$ 和 ${\mu }_w$ 平均值时的剪切速率， $P_a$ 为经济系数。

8) 残余饱和度

各相残余饱和度与毛管数 $N_c$ 有关，通过实验可测得不同毛管数下的残余饱和度。

$N_c=\frac{\left|{\displaystyle \sum {}_j{\mu }_j{\bar{V}}_j}\right|}{\sigma },$ (60)

$S_{rj}=S_{rj}\left(N_c\right),$ (61)

式中 $S_{rj}$ 表示j相的残余饱和度。