引言

目前中国各大油田的主力油区都已进入高含水开发阶段，化学驱技术已成为提高油田采收率的重要方法并开始大规模应用。聚合物驱油是目前技术最为成熟，并已大规模工业化应用的一种化学驱技术。在表面活性剂-碱-聚合物三元复合驱油体系中，聚合物同样是重要组分，主要起到改善流度比的作用^[1]。同水驱相比，聚合物驱油有很多新的物理化学现象，如黏度与浓度关系、剪切变稀、聚合物吸附以及渗透率下降等，聚合物水溶液在地下运移机理也有很大区别^[2]。 为正确评价聚合物驱后地层参数、聚合物堵水效果，现场往往进行大量的压力测试，并通过对测试资料的解释，判断井的特性参数，污染程度，推算地层参数，为调整开发方案服务。因此，为满足当前试井分析需要，多相、多组分的聚合物驱数值试井模型的开发势在必行。

目前聚合物驱试井方法主要有解析试井和数值试井^[3-5]。在解析试井方面，Ikoku C U^[5]和Odeh A S等^[6]分别较早地给出了非牛顿幂律流体在均质多孔介质中流动的数学模型，并给出了解析解。Vongvuthipornchai S和Raghavan R^[7]给出了考虑井储和表皮效应的非牛顿流试井典型曲线。中国，宋考平等^[8]采用复合油藏假设对非牛顿幂律流体驱油的渗流规律进行了理论研究，给出了井底压力及导数曲线，张传宝等^[9] 研究了非牛顿幂律流体的试井解释方法。由于解析试井在处理复杂边界、多相流动等复杂问题上的局限性以及忽略了聚合物驱中复杂的物理化学现象，并且只能将聚合物溶液简化为幂律流体等原因，数值试井是近年来主要发展方向。

Mahani H等^[10]发展了一种解析解与数值解结合的非牛顿流试井解释方法，并通过压力降落测试数据对聚合物非牛顿特性进行反演。在中国，吴明录等^[11]根据流线模拟方法建立了聚合物驱流线数值试井解释模型。总的来说，国内外关于非牛顿流体数值试井方面的研究较少，目前还缺乏能对实际聚合物驱油藏进行试井分析解释的实用模拟器和成熟的试井解释方法。

国外基于PEBI网格的数值试井已形成商业化软件，但只是基于黑油模型的三相流模块，没有聚合物驱数值试井模块。商业化的数值模拟软件包括聚合物驱模块，但数值模拟软件计算时间步长大，不能对开关井瞬态压力进行模拟，也没有考虑井筒存储，所以聚合物驱数值模拟软件不能计算瞬时井底压力。

本文基于PEBI网格^[12-14]建立了两相(油、水)四组分(油、水、聚合物、阳离子)的全隐式聚合物驱数值试井模型，考虑了聚合物流变性、吸附、不可及孔隙体积、渗透率下降系数及残余阻力系数等聚合物驱物理化学现象和参数。研究了剪切变稀效应对试井曲线的影响，并用实测数据进行了验证。

1 物化现象

与水驱相比，聚合物驱过程中会产生一系列新的物理化学现象，主要包括非牛顿性、渗透率下降、聚合物吸附、扩散弥散、不可及孔隙体积以及阳离子对聚合物溶液的影响等。

注入地下的聚合物可以在岩石表面发生吸附，也可在油藏孔隙介质中产生滞留，从而引起聚合物的消耗损失。由于高相对分子质量聚合物的体积较大，使得它难以进入油藏中一些较小的孔隙，这对聚合物来说相当于有效孔隙度的降低。聚合物溶液在通过多孔介质时，除了由黏度产生的流度变化外，由于聚合物的吸附滞留，也会引起流度的减小，用渗透率下降系数来表示，在聚合物溶液通过多孔介质后，由于聚合物的吸附滞留，对聚合物后的水驱产生残余阻力，使得水相流度比聚合物溶液通过之前要小，通过残余阻力系数来描述^[15]。

如果油藏为等温模型，则聚合物溶液的黏度与溶液中包含的聚合物、阳离子等各个组分浓度相关。黏度的计算，可以使用现有的公式，即给出聚合物黏度与聚合物浓度、盐浓度之间的近似解析关系式，也可以通过给出基于实验的数据表进行插值计算。

聚合物溶液是非牛顿流体，黏度会随剪切率变化而改变。剪切作用越强，黏度越小。一定浓度下聚合物溶液黏度与剪切率之间的关系可以用公式(1)来表示(Meter方程)

$\mu = \mu _\infty + \dfrac{{\mu _{0} - \mu _\infty }}{{1 + \left( {\dfrac{{\dot \gamma }}{{\dot \gamma _{1/2} }}} \right)^{P_\alpha - 1} }}$

(1)

式中：$\mu _\infty$—剪切率趋于无穷大时的聚合物溶液黏度，mPa$\cdot$s；

$\mu $—聚合物流动时的有效黏度，mPa$\cdot$s；

$\mu _{0}$—零剪切速率下的黏度，mPa$\cdot$s；

$P_\alpha$—无因次常数；

$\dot \gamma$—剪切率，s$^{-1}$；

$\dot \gamma _{1/2}$—$\mu = (\mu _{\rm{o}} - \mu _{\rm{w}} )/2$时所对应的剪切率，s$^{-1}$。

$P_\alpha$，$\dot \gamma$等参数难以确定，因而常用数据表形式给出黏度随速度的变化关系。可以用公式(2)来描述剪切作用下的聚合物溶液黏度。

$\mu = M\mu _0 + (1 - M)\mu _{\rm{w}}$

(2)

式中： $\mu _{\rm{w}} $—水的黏度，mPa$\cdot$s；

M—随渗流速度变化的剪切变稀系数，无因次，介于0$\sim$1。

2 数学模型

中国老油田地下主要为油水两相流动，采用聚合物驱开采时，水相黏度主要受聚合物和盐浓度影响。因此，本文采用的聚合物驱数学模型为两相(油、水)，四组分(油、水、聚合物、一价阳离子)。

聚合物驱数学模型的基本假设条件为：(1) 油藏为等温模型；(2) 扩展的达西定律适合于描述多相流动；(3) 流体由油、水两相和油、淡水、聚合物、一价阳离子4个组分组成；(4) 聚合物和一价阳离子存在于水相中；(5) 聚合物在岩石表面的吸附-滞留为不可逆过程；(6) 除了聚合物的吸附外没有其他化学反应的发生。

根据上述假设，聚合物驱的基本方程如下

油

$\nabla \left[{\dfrac{{KK_{{\rm{ro}}} }}{{\mu _{\rm{o}} B_{\rm{o}} }}\left( {\nabla p_{\rm{o}} - \gamma _{\rm{o}} \nabla D} \right)} \right] + q_{\rm{o}} = \dfrac{\partial }{{\partial t}}\left( {\dfrac{{\phi S_{\rm{o}} }}{{B_{\rm{o}} }}} \right)$

(3)

水

$\nabla \left[{\dfrac{{KK_{{\rm{rw}}} }}{{R_{\rm{k}} \mu _{\rm{w,eff}} B_{\rm{w}} }}\left( {\nabla p_{\rm{w}}\!-\!\gamma _{\rm{w}} \nabla D} \right)} \right]\!+\!q_{\rm{w}}\!=\! \dfrac{\partial }{{\partial t}}\left( {\dfrac{{\phi S_{\rm{w}} }}{{B_{\rm{w}} }}} \right)$

(4)

聚合物

$\nabla \left[{\dfrac{{KK_{{\rm{rw}}} }}{{R_{\rm{k}} \mu _{\rm{w,eff}} B_{\rm{w}} }}\left( {\nabla p_{\rm{w}} - \gamma _{\rm{w}} \nabla D} \right)C_{\rm{p}} } \right] + q_{\rm{p}} = \\ \dfrac{\partial }{{\partial t}}\left[{\dfrac{{\phi _{\rm{p}} \left( {S_{\rm{w}} C_{\rm{p}} + C_{{\rm{pads}}} } \right)}}{{B_{\rm{w}} }}} \right]$

(5)

一价阳离子

$\nabla \left[{\dfrac{{KK_{{\rm{rw}}} }}{{R_{\rm{k}} \mu _{\rm{w,eff}} B_{\rm{w}} }}\left( {\nabla p_{\rm{w}} - \gamma _{\rm{w}} \nabla D} \right)C_{{\rm{b}}} } \right] + q_{{\rm{cl}}} = \\ \dfrac{\partial }{{\partial t}}\left( {\dfrac{{\phi S_{\rm{w}} C_{{\rm{b}}} }}{{B_{\rm{w}} }}} \right)$

(6)

饱和度关系

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

(7)

毛管压力方程

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

(8)

井方程

$q_{{\rm{lsc}}} = \dfrac{1}{{\mu _l B_l }}\dfrac{{\theta K_{{\rm{r}}l} Kh}}{{\ln \left( {r_{\rm{e}} /r_{\rm{w}} } \right) + S}}\left( {p - p_{{\rm{wf}}} } \right)$

(9)

式中： K—绝对渗透率，D；

$K_{\rm{r}}$—相对渗透率，无因次；

B—地层体积系数；

$\mu _{\rm{w,eff}}$—水相有效黏度，Pa$\cdot$s；

$\phi$—孔隙度，无因次；

t—时间，s；

q—产量，m³/d；

$S_l$—相饱和度，l=o，w；

p—压力，Pa；

$\gamma$—重度，N/m³；

D—深度，m；

$R_{\rm{k}}$—渗透率下降系数，无因次；

$C_{\rm{p}}$—聚合物浓度，kg/m³；

$C_{{\rm{pads}}} $—聚合物吸附浓度，kg/m³；

$C_{\rm{cl}}$—一价阳离子浓度，kg/m³；

$\sigma$—表面张力系数，N/m；

$p_{\rm{c}}$—毛管压力，Pa；

$p_{\rm{cow}}$—油水界面毛管压力，Pa；

$q_{\rm{lsc}}$— l相地面产量，m³/d；

$p_{\rm{wf}}$—井底流压，Pa；

S—表皮因子，无因次；

$\theta$—井筒射开角度，rad；

$r_{\rm{w}}$—井筒半径，m；

$r_{\rm{e}}$—等效半径，m；

h—地层厚度，m；

下标

w，o，p，cl—水，油，聚合物，一价阳离子。

3 数值模型

由于求解区域较大，而井半径仅0.1 m。同时地层往往存在断层及复杂边界等。基于PEBI 网格采用有限体积法离散，并形成一套完整的网格划分及数值求解模型^[16]。以聚合物方程为例，给出其离散格式

$\begin{align} & \sum\limits_{j}{{{\left[ {{T}_{ij,\text{w}}}{{C}_{\text{p}}}\left( \Delta p-{{\gamma }_{\text{w}}}\Delta D \right) \right]}^{n+1}}}={{C}_{\text{pp}}}\delta p+{{C}_{\text{pw}}}\delta {{S}_{\text{w}}}+ \\ & {{C}_{\text{pc}}}\delta {{C}_{\text{p}}}+{{C}_{\text{pa}}}\delta {{C}_{\text{ads}}}+q_{\text{wsc}}^{n+1} \\ \end{align}$

(10)

式中：

$T_{ij,{\rm{w}}} = \left( {\dfrac{{KK_{{\rm{rw}}} }}{{\mu _{{\rm{w,eff}}} B_{\rm{w}} R_{\rm{k}} }}} \right)_{ij} \dfrac{{w_{ij} }}{{d_{ij} }}{；}$

$C_{{\rm{pp}}} = \dfrac{{V_i }}{{\Delta t}}\left[{\dfrac{1}{{B_{\rm{w}}^n }}\dfrac{{\partial \phi _{\rm{p}} }}{{\partial p}}\left( {C_{\rm{p}}^n S_{\rm{w}}^n + C_{{\rm{pads}}}^n } \right) + \phi _{\rm{p}}^{n + 1} \dfrac{{\partial \left( {1/B_{\rm{w}} } \right)}}{{\partial p}}\left( {C_{\rm{p}} S_{\rm{w}} + C_{{\rm{pads}}} } \right)} \right]{；}$

$C_{{\rm{pw}}} = \dfrac{{V_i }}{{\Delta t}}\left( {\dfrac{{\phi _{\rm{p}} C_{\rm{p}} }}{{B_{\rm{w}} }}} \right)^{n + 1}{；}$

$C_{{\rm{pc}}} = \dfrac{{V_i }}{{\Delta t}}\left( {\dfrac{{\phi _{\rm{p}} }}{{B_{\rm{w}} }}} \right)^{n + 1} S_{\rm{w}}^n{；}$

$C_{{\rm{pa}}} = \dfrac{{V_i }}{{\Delta t}}\left( {\dfrac{{\phi _{\rm{p}} }}{{B_{\rm{w}} }}} \right)^{n + 1}{；}$

$\Delta D$—相邻网格点纵向坐标差，$\Delta D$=$D_j-D_i$；

w—网格交界面面积，m²；

d—相邻网格中心点之间的距离，m；

V—网格体积，m³；

n—时间步；

下标

i—当前网格编号；

j—相邻网格编号。

全隐式线性化

$\begin{align} & \sum\limits_{j}{\left[ T_{ij.\text{w}}^{\nu }\left( \delta {{p}_{j}}-\delta {{p}_{i}}+p_{j}^{\nu }-p_{i}^{\nu }-{{\gamma }_{\text{w}}}\Delta D \right)+ \right.} \\ & {{\left( \frac{\partial {{T}_{ij.\text{w}}}{{C}_{\text{p}}}}{\partial p} \right)}^{v}}+\delta p+\left( p_{j}^{\nu }-p_{i}^{\nu }-{{\gamma }_{\text{w}}}\Delta D \right) \\ & \left. +{{\left( \frac{\partial {{T}_{ij.\text{w}}}{{C}_{\text{p}}}}{\partial {{C}_{\text{p}}}} \right)}^{v}}+\delta {{C}_{p}}+\left( p_{j}^{\nu }-p_{i}^{\nu }-{{\gamma }_{\text{w}}}\Delta D \right) \right] \\ & =C_{\text{pp}}^{\nu +1}\left( \delta p+{{p}^{\nu }}-{{p}^{n}} \right)+C_{\text{pw}}^{\nu +1}\left( \delta {{S}_{\text{w}}}+{{S}_{w}}-S_{\text{w}}^{n} \right)+ \\ & {{C}_{\text{pc}}}\left( \delta {{C}_{\text{p}}}+C_{\text{p}}^{\nu }-C_{\text{p}}^{n} \right)+{{C}_{\text{pa}}}\left( \delta {{C}_{\text{pads}}}+C_{\text{pads}}^{\nu }-C_{\text{pads}}^{n} \right)+q_{\text{wsc}}^{n+1} \\ \end{align}$

(11)

同样以聚合物组分为例，源汇相离散方程为

$\begin{align} & {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu +1}}={{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{p}_{\text{w},i}}}\delta {{p}_{i}} \\ & +\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{S}_{\text{w}}}}\delta {{S}_{\text{w}}}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{C}_{\text{p}}}}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{C}_{cl}}}\partial {{C}_{cl}}+ \\ & \frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{p}_{wfref,w}}}\partial {{p}_{wfref,w}} \\ \end{align}$

(12)

对定流量生产有

$\begin{align} & \frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{p}_{\text{w},i}}}\delta {{p}_{i}}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{S}_{\text{w}}}}\delta {{S}_{\text{w}}}+ \\ & \frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{C}_{\text{p}}}}\partial {{C}_{\text{p}}}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{C}_{cl}}}\partial {{C}_{cl}}+\frac{\partial {{\left( {{q}_{\text{ws},i}}{{C}_{\text{p}}} \right)}^{\nu }}}{\partial {{p}_{wfref,w}}}\partial {{p}_{wfref,w}}=0 \\ \end{align}$

(13)

式中： $\nu$—迭代步；

$\delta f={{f}^{\begin{matrix} (v+1) \\ n+1 \\ \end{matrix}}}-{{f}^{\begin{matrix} (v) \\ n+1 \\ \end{matrix}}}={{f}^{v+1}}-{{f}^{v}}-f,f={{p}_{j}},{{p}_{i}},{{C}_{p}},{{S}_{w}},{{C}_{pads}};$ $q_{\text{wsc}}^{n+1}$—水相在标准条件下的产量，m³/d； $p_{{\rm{wfref,w}}}$—参考井底压力，Pa；

下标

+—上游网格的属性。

4 计算结果及分析 4.1 正确性验证

为了验证本模型的正确性，对一注一采的注入井井底压力进行数值模拟，并与目前主流商业油藏数值模拟软件Eclipse计算结果进行对比。

本文的计算模型在将井视作源或汇，网格变成矩形网格({图 1a) 时可以进行油藏数值模拟，其模拟结果与商用软件一致。计算区域为200 m$\times$200 m 的封闭油藏区域。油藏中有两口井(图 1a)，井1为注入井，井2为生产井，表 1给出了主要输入参数，图 2给出了相渗曲线。

图 3是聚合物浓度为0时注入井井网格压力计算结果与Eclipse对比，图 4为注入聚合物浓度$C_{\rm{p}}$=0.005，0.010 kg/m³时的注入井井网格压力计算结果与Eclipse的对比。

图 3表明，在聚合物浓度为零、聚驱模型退化为黑油模型时，模型对井底压力计算结果与Eclipse计算结果基本一致。图 4表明，在聚合物浓度不为零时，二者计算结果基本平行，压力值相对误差很小。对比结果验证了本文模型的正确性。

4.2 剪切变稀曲线对井底瞬态压力的影响

剪切变稀是聚合物驱中重要流体性质，本文采用3种不同的剪切变稀曲线进行数值模拟。图 5给出了3条剪切变稀曲线，其中曲线1表示的聚合物溶液剪切变稀性质最弱，曲线3表示剪切变稀作用最强，曲线2所表示的剪切变稀特性介于曲线1和3之间。

仍采用4.1中两井(一注一采)算例进行计算。改变部分参数：水平渗透率为0.5 mD，孔隙度为0.36，并改变生产制度：生产井以70 m³/d的产量进行生产。注入井以70 m³/d的注入量注入聚合物溶液，关井110 d。为正确模拟井底瞬时压力，采用可变PEBI网格(图 1b)计算，在井处理上采用peaceman模型，以提高计算速度。图 6给出了剪切变稀曲线1，2，3下和无剪切变稀作用时计算得到的注入井关井后井底压力变化及其导数曲线。

从图 6中看到：考虑剪切变稀后，井底压力降落及导数曲线有明显的变化。首先，考虑剪切变稀后，注入井压力降落变小。这是因为剪切作用越大，黏度越小，流动能力越强，所以压力降落和压力导数越小。其次，考虑剪切变稀后，压力导数前期降低后期重合，并且在井储段结束后，压力导数出现“V”形抖动。出现这种现象的原因是，关井后较短时间内注入井周围流动速度急剧降低，由于聚合物的非牛顿剪切变稀特性，水相黏度会大幅增加，流动能力突然变弱，压力恢复变慢，导数曲线出现下凹。图 7为注入井关井后井周围水相速度变化，从图 7中可以看到，在无剪切作用以及剪切变稀作用很小时，速度降低较为平缓；而在剪切变稀作用较强时，关井后水相速度变化则非常急剧，并且还注意到图 7中速度曲线急剧变化的时间区间与图 6中导数曲线出现抖动的时间区间是相同的，印证了关井后井周围水相速度的急剧降低和聚合物剪切变稀特性是导致导数曲线抖动的原因。

还通过数值计算研究了在考虑剪切变稀的情况下主要地层参数(渗透率)对井底瞬态压力的影响。图 8 给出了剪切曲线3下，改变地层渗透率计算得到的井底压力及导数曲线。从图中可以看到，随着渗透率增大，出现下凹的时间段向左下平移。这是因为渗透率越大，速度下降越快，聚合物黏度变化快，导致下凹的时间段前移。图 8显示渗透率变化对压力导数曲线抖动的幅度影响不大。

4.3 考虑剪切变稀作用的试井拟合方法

由上可知，聚合物溶液在地层中的剪切变稀特性和地层渗透率对聚驱试井瞬态压力曲线影响很大。聚合物溶液的剪切变稀作用使导数曲线在早期阶段出现有规律的抖动。地层渗透率对瞬态压力曲线的影响与水驱规律基本相同，在聚驱情况下还会导致导数曲线抖动的位置发生移动，而对抖动的幅度影响很小。本文结果对于聚驱瞬态压力测试数据分析具有重要指导意义，给出了考虑剪切变稀的聚合物驱数值试井拟合的几个最主要步骤：

(1) 在双对数图上作出实测压力降落和压力导数曲线。

(2) 通过后期曲线形态分析调整油藏边界及地层渗透率。

(3) 通过早期阶段导数曲线的拟合判断剪切变稀效应，由于压力导数在此阶段对剪切变稀曲线的变化非常敏感，可通过对导数曲线抖动形状的拟合提高拟合精度，关井后水相速度变化曲线有助于剪切变稀曲线的合理调整。

(4) 拟合井储段曲线确定井储。

(5) 数值试井拟合过程为“试错法”，即反复执行步骤(2)$\sim$(4)直到满足拟合精度要求。

5 实例分析

为了验证考虑剪切变稀作用对于解释聚驱现场压力测试资料的实际意义，在大庆油田注聚合物区块一口实际油井中关井并下压力计至井底，测量12 h的井底压力随时间变化曲线。选用一个五点井网进行模式计算，4口注入井以60 m³/d的总注入量注入浓度为1.2 kg/m³的聚合物溶液。生产井产量为60 m³/d，生产84.5 d。采用油田实测的相渗曲线，聚合物的相关参数采用油田提供的实验数据。

图 9给出了理论计算与实际测试的井底压力及其导数曲线图，从图中可以看出，计算的压力降落与实测数据基本重合。实测与计算的压力导数曲线形态一致，基本拟合，说明剪切变稀是聚驱压力导数“V”形抖动的内在原因。

6 结 论

（1） 建立了两相、四组分的聚合驱数学模型，考虑了聚合物溶液黏度、剪切变稀、聚合物吸附、残余阻力系数、不可及孔隙体积等物理化学现象，并基于PEBI网格进行了数值求解，研发了聚合物驱数值试井计算模型

（2） 研究地层和流体参数对聚驱关井后瞬态压力曲线动态响应的影响，发现聚驱井底压力导数曲线在早期阶段会出现明显的抖动，结合关井后水相速度变化分析得出了聚合物的剪切变稀特性是导致压力导数曲线早期抖动的结论。数值计算的结果表明，地层渗透率的变化会导致导数曲线抖动的位置发生明显变化但对抖动的幅度影响很小。基于以上结论给出了考虑剪切变稀作用的聚驱数值试井拟合方法。

（3） 对大庆油田某注聚井实测压力数据的拟合，验证了数值模拟所发现规律的正确性，其对现场压力测试资料解释、地层参数反演有很大指导意义。