目前，水驱开发是油气藏主要开发方式，对水驱油藏来说，油井生产历史准确拟合和产量准确预测都有赖于油藏油水相对渗透率关系和Leverett分流函数^[1-2]。目前，人们广泛应用相对渗透率比值和流体饱和度函数分析油气藏生产动态和产水规律，揭示油气藏中油水渗流关系^[3-7]。并且，主要靠实验室一维油水驱替实验，结合Darcy方程，根据JBN方法或其改进方法，得到相对渗透率与流体饱和度函数关系^[8-14]。然而，实际油藏存在较强非均质性，实验室一维岩心与油藏存在较大差距，并不能代表油藏非均质情况；同时，水压驱动下，油水井近井地带为径向流，油藏深处方为单向渗流。因此，实验室一维单向渗流并不能真实地反映近井地带油水渗流关系。计秉玉运用油水相对渗透率曲线和径向Bockley-Leverett方程确定油水总流度径向变化的方法，建立了油水两相径向稳定渗流条件下压力分布和产量计算公式^[15]，但是他应用的相对渗透率曲线仍然是一维油水驱替实验得到的。秦积舜等利用平面物理模型研究了油水两相平面渗流的拟相渗特征^[16]，但是计算方法仍然是常规JBN方法。为此，基于平面径向流动模型，笔者建立了平面径向流渗流方程，求得了油水相对渗透率比值与含水饱和度解析解，并用上述思路分析了平面径向流实验结果，用相对渗透率比值描述两相径向流中油相和水相渗流关系。与一维渗流模型相比，平面径向流相对渗透率比值与流体饱和度呈直线关系的饱和度区间变小。渗流中，油相和水相之间的相互作用以非线性作用为主，这无疑对认识油气藏近井地带油水两相径向渗流特征，准确模拟和预测油水井动态，提供了新观点和新证据。

1 油水径向流渗流理论 1.1 平面径向渗流方程建立

砂岩模型如图 1所示。砂岩模型均质，孔隙度φ、水相渗透率k_w、模型厚度h，中心井井径r_w、模型半径r_e，模拟油水等流体向中心井渗流的过程。

假设油水不可压缩，非混相；并且忽略毛管力对渗流过程影响，即p_c=0，水注入速度为q_i，边缘流体压力为p_e，井底流压p_w，流动压差则为Δp=p_e-p_w。

水驱油过程，注水速度与产液速度相等。

根据Darcy方程，中心井产油和产水速度为

$ {q_{\rm{w}}} = - \frac{{2{\rm{ \mathsf{ π} }}hrk{k_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}}\frac{{{\rm{d}}p}}{{{\rm{d}}r}}, $

(1)

$ {q_{\rm{o}}} = - \frac{{2{\rm{ \mathsf{ π} }}hrk{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}}}\frac{{{\rm{d}}p}}{{{\rm{d}}r}}。$

(2)

流体流量为

$ {q_{\rm{i}}} = - 2{\rm{ \mathsf{ π} }}hrk\left( {\frac{{{k_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}} + \frac{{{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}}}} \right)\frac{{{\rm{d}}p}}{{{\rm{d}}r}}。$

(3)

定义总相对流度λ_r，

$ {\lambda _{\rm{r}}} = \frac{{{k_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}} + \frac{{{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}}} = {\lambda _{{\rm{rw}}}} + {\lambda _{{\rm{ro}}}}。$

(4)

那么

$ - {\rm{d}}p = \frac{{{q_{\rm{i}}}}}{{2{\rm{ \mathsf{ π} }}hk{\lambda _{\rm{r}}}}}\frac{{{\rm{d}}r}}{r}。$

(5)

积分，得到压差与时间的关系：

$ \Delta p = {p_{\rm{e}}} - {p_{\rm{w}}} = - \frac{{{q_{\rm{i}}}}}{{2{\rm{ \mathsf{ π} }}hk}}\int_{{r_{\rm{w}}}}^{{r_{\rm{e}}}} {\frac{{{\rm{d}}r}}{{r{\lambda _{\rm{r}}}}}} 。$

(6)

方程边界条件为

$ p\left| {_{r = {r_{\rm{e}}}}} \right. = {p_{\rm{e}}},p\left| {_{r = {r_{\rm{w}}}}} \right. = {p_{\rm{w}}}。$

(7)

方程(6)、(7)即为两相径向渗流数学模型。

1.2 两相径向流方程求解

总相对流度λ_ro为饱和度S_w的函数，将式(6)做积分变换，成为S_w函数，边界条件变为

$ \begin{array}{*{20}{c}} {{f_{\rm{w}}}\left| {_{r = {r_{\rm{e}}}}} \right. = 1,{{f'}_{\rm{w}}}\left| {_{r = {r_{\rm{e}}}}} \right. = 0,{\rm{d}}{{f'}_{\rm{w}}}\left| {_{r = {r_{\rm{e}}}}} \right. = 0,}\\ {{f_{\rm{w}}}\left| {_{r = {r_{\rm{w}}}}} \right. = {f_{{\rm{w}}2}},{{f'}_{\rm{w}}}\left| {_{r = {r_{\rm{w}}}}} \right. = }\\ {{{f'}_{{\rm{w2}}}}\left| {_{r = {r_{\rm{w}}}}} \right. = - \frac{{2{\rm{ \mathsf{ π} }}h\phi r_{\rm{w}}^2}}{{{V_{\rm{i}}}\left( t \right)}},{\rm{d}}{{f'}_{\rm{w}}}\left| {_{r = {r_{\rm{w}}}}} \right. = }\\ {{\rm{d}}{{f'}_{\rm{2}}}\left| {_{r = {r_{\rm{w}}}}} \right. = \frac{{2{\rm{ \mathsf{ π} }}h\phi r_{\rm{w}}^2}}{{V_{\rm{i}}^2\left( t \right)}}{q_{\rm{i}}}{\rm{d}}t。} \end{array} $

(8)

径向渗流情形下，Buckley-Leverrt驱替方程的通解为

$ {r^2}\left( {{S_{\rm{w}}}} \right) = - \frac{{{V_{\rm{i}}}\left( t \right)}}{{2{\rm{ \mathsf{ π} }}h\phi }}f'\left( {{S_{\rm{w}}}} \right), $

(9)

其中，$f'\left( {{S_{\rm{w}}}} \right) = \frac{{{\rm{d}}{f_{\rm{w}}}}}{{{\rm{d}}{S_{\rm{w}}}}};{V_{\rm{i}}}\left( t \right) = \int_0^t {{q_i}{\rm{d}}t} $。

将式(9)微分，得到

$ \frac{{{\rm{d}}r}}{r} = \frac{1}{2}\frac{{{\rm{d}}{{f'}_{\rm{w}}}}}{{{{f'}_{\rm{w}}}}}。$

(10)

将式(10)带入到式(6)，并根据边界条件(8)，变换积分限，得$\Delta p = - \frac{{{q_i}}}{{4\pi kh}}\int_{{f_{{\rm{w}}2}}}^0 {\frac{{{\rm{d}}{f_{\rm{w}}}'}}{{{f_{\rm{w}}}'{\lambda _{\rm{r}}}}}} $。

将上述方程引入流函数G，即

$ \begin{array}{l} \Delta p = - \frac{{{q_{\rm{i}}}}}{{4{\rm{ \mathsf{ π} }}kh}}\int_{{f_{{\rm{w2}}}}}^0 {\frac{{{\rm{d}}{{f'}_{\rm{w}}}}}{{{{f'}_{\rm{w}}}{\lambda _{\rm{w}}}}} = } \\ \;\;\;\;\frac{{{q_{\rm{i}}}}}{{4{\rm{ \mathsf{ π} }}kh}}\int_0^{{f_{{\rm{w2}}}}} {\frac{{{\rm{d}}{{f'}_{\rm{w}}}}}{{{{f'}_{\rm{w}}}{\lambda _{\rm{r}}}}}} = \frac{{{q_{\rm{i}}}}}{{4{\rm{ \mathsf{ π} }}kh}}G。\end{array} $

(11)

根据式(11)，可以得到$G = - \frac{{4\pi kh\Delta p}}{{{q_i}}} = \int_0^{{f_{{\rm{w}}2}}} {\frac{{{\rm{d}}{f_{\rm{w}}}'}}{{{f_{\rm{w}}}'{\lambda _{\rm{r}}}}}} $。这里，流函数G是r、f′_w的函数，与Buckley-Leverrt线性驱替中流函数不同，不仅依赖于含水率变化，还受径向渗流面积的影响，该参数揭示了径向流和Buckley-Leverrt流之间的差异，线性驱替中，该函数与r无关。

将式(11)对时间微分

$ \frac{{{\rm{d}}\left( {\Delta p} \right)}}{{{\rm{d}}t}} = \frac{1}{{4{\rm{ \mathsf{ π} }}kh}}\left[ {{q_{\rm{i}}}\frac{{{\rm{d}}G}}{{{\rm{d}}t}} + G\frac{{{\rm{d}}{q_{\rm{i}}}}}{{{\rm{d}}t}}} \right], $

(12)

式中：${q_{\rm{i}}} = \frac{{{\rm{d}}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}t}},\frac{{{\rm{d}}{q_{\rm{i}}}}}{{{\rm{d}}t}} = \frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}}$。

所以

$ G\frac{{{\rm{d}}{q_{\rm{i}}}}}{{{\rm{d}}t}} = \frac{{4{\rm{ \mathsf{ π} }}hk\Delta p}}{{{q_{\rm{i}}}}}\frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}}。$

(13)

根据边界条件，得到

$ {{f'}_{{S_{{\rm{w2}}}}}} = \frac{{q_{\rm{i}}^2}}{{{\lambda _{\rm{r}}}\left( {{S_{{\rm{w2}}}}} \right){V_{\rm{i}}}\left( t \right)}}。$

(14)

那么流动压差随时间变化特征为

$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}\Delta p}}{{{\rm{d}}t}} = \frac{1}{{4{\rm{ \mathsf{ π} }}kh}}\left[ {\frac{{4{\rm{ \mathsf{ π} }}hk\Delta p}}{{{q_{\rm{i}}}}}\frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}} - } \right.}\\ {\left. {\frac{{q_{\rm{i}}^2}}{{{\lambda _{\rm{r}}}\left( {{S_{{\rm{w2}}}}} \right){V_{\rm{i}}}\left( t \right)}}} \right] = \frac{{\Delta p}}{{{q_{\rm{i}}}}}\frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}} - }\\ {\frac{{q_{\rm{i}}^2}}{{4{\rm{ \mathsf{ π} }}kh{\lambda _{\rm{r}}}\left( {{S_{{\rm{w2}}}}} \right){V_{\rm{i}}}\left( t \right)}}。} \end{array} $

(15)

根据物理模拟实验条件，将恒压和恒速两种边界条件应用于流动方程。

a) 恒压条件下，Δp为常数，V_i(t)为变量，那么式(15)左边为0，总相对流度λ_r则为

$ {\lambda _{\rm{r}}}\left( {{S_{{\rm{w2}}}}} \right) = \frac{{q_{\rm{i}}^3}}{{4{\rm{ \mathsf{ π} }}kh{V_{\rm{i}}}\left( t \right)\Delta p\frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}}}}。$

(16a)

b) 恒速条件下，q_i为常数，Δp为变量，此时，

$ \frac{{{{\rm{d}}^2}{V_{\rm{i}}}\left( t \right)}}{{{\rm{d}}{t^2}}} = 0,{\lambda _{\rm{r}}}\left( {{S_{{\rm{w2}}}}} \right) = \frac{{{q_{\rm{i}}}}}{{4{\rm{ \mathsf{ π} }}kh{V_{\rm{i}}}\left( t \right)\frac{{{\rm{d}}\Delta p}}{{{\rm{d}}t}}}}。$

(16b)

1.3 相对渗透率计算

根据分流量方程，${f_{\rm{w}}} = \frac{{{q_{\rm{w}}}}}{{{q_{\rm{i}}}}} = \frac{{{k_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}{\lambda _{\rm{r}}}}}$、${f_{\rm{o}}} = \frac{{{q_{\rm{o}}}}}{{{q_{\rm{i}}}}} = \frac{{{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}{\lambda _{\rm{r}}}}}$得

$ \frac{{{k_{{\rm{rw}}}}}}{{{k_{{\rm{ro}}}}}} = \frac{{{f_{\rm{w}}}{\mu _{\rm{w}}}}}{{{f_{\rm{o}}}{\mu _{\rm{o}}}}}。$

(17)

因此，根据式(17)和式(16a)或(16b)即可得到相对渗透率。

1.4 饱和度方程计算

实验中，S_wi为束缚水饱和度，当注水量为V_i(t)时，产油量为V_o(t)，产水量V_w(t)，因此

$ {V_{\rm{i}}}\left( t \right) = {V_{\rm{w}}}\left( t \right) + {V_{\rm{o}}}\left( t \right)。$

(18)

并且，产量递减模型^[4]

$ \frac{{{V_{\rm{i}}}}}{{{V_{\rm{o}}}}} = a + b\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}。$

(19)

水相在径向岩心介质中平均饱和度为

$ {{\bar S}_{\rm{w}}} = {S_{{\rm{wi}}}} + \frac{{{V_{\rm{o}}}}}{{{V_{\rm{p}}}}}。$

(20)

根据Toth等研究^[8]，水相在界面上饱和度为

$ {S_{{\rm{w2}}}} = {S_{{\rm{wi}}}} + b{\left[ {\frac{{\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}}}{{a + b\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}}}} \right]^2}。$

(21)

含油率、含水率和径向流水驱效率为

$ {f_{\rm{o}}} = \frac{a}{{{{\left[ {a + b\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}} \right]}^2}}},{f_{\rm{w}}} = 1 - {f_{\rm{o}}}, $

(22)

$ {E_{\rm{D}}} = \frac{{{V_{\rm{o}}}}}{{{V_{\rm{p}}}\left( {1 - {S_{{\rm{wi}}}}} \right)}} = \frac{{\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}}}{{\left( {a + b\frac{{{V_{\rm{i}}}}}{{{V_{\rm{p}}}}}} \right)\left( {1 - {S_{{\rm{wi}}}}} \right)}}。$

(23)

2 室内平面径向渗流实验及分析 2.1 岩心物性参数及实验方法

实验用径向砂岩模型如图 1所示，岩心参数和流体参数见表 1。

实验过程如下：首先将岩心抽真空，饱和地层水。然后将岩心装入驱替流程，如图 2所示。岩心顶面和底面密封，侧面与空气相通。注入流体从岩心顶部的中心孔注入，被驱出的流体从岩心外侧面流出，然后通过装置底部的管线流入量筒计量。压差通过压力传感器测定。首先从中间的入口注入油相进行油驱水，直至出口端不再含水，模拟油藏形成过程，建立束缚水或初始含水饱和度；然后从中间的入口注入水相，模拟油藏水驱开发过程。此时，径向模型中油水渗流特征则反映油藏近井地带油水渗流特点。

实验中，采用恒压驱替，水驱油径向压差为0.33 MPa，水沿边缘注入，实验中记录注水量、中心井产油量、产水量。图 3为根据方程(19)得到实验结果。

计算参数a、b，根据式(21)计算水相饱和度，并结合式(22)，得到f_w-S_w关系，同时，根据Buckley-Leverrt一维渗流模型，计算相应的f_w-S_w关系曲线，计算结果与实验结果对比见图 4。图 4表明，与Buckley-Leverrt一维渗流模型相比，笔者得到的径向流渗流模型与实际实验中的油水渗流过程更为接近，更能反映径向渗流过程中油水渗流关系。

根据式(23)可以计算径向渗流水驱油效率，它与基于Buckley-Leverrt一维渗流模型得到的驱油效率、实验结果对比见图 5。

图 5表明，与Buckley-Leverrt一维渗流模型相比，这里得到的径向流渗流模型与实际径向流驱油结果更为接近。

2.2 径向流模型与一维线性流模型计算结果对比分析

根据式(17)、(21)和(22)计算径向渗流过程中油相和水相相对渗透率比值，为了对比，还利用传统的JBN方法对数据进行了计算，结果如图 6所示。

在水驱油藏中，油相和水相在各自的流动通道系统中渗流，相对渗透率比值主要是含水饱和度线性函数(在单对数坐标上)，对应着两相流范围。

描述方程为：

$ \frac{{{k_{{\rm{ro}}}}}}{{{k_{{\rm{rw}}}}}} = \frac{{{k_{\rm{o}}}}}{{{k_{\rm{w}}}}} = A{{\rm{e}}^{ - B{S_{\rm{w}}}}}。$

(24)

式中，系数A、B为相对渗透率比值曲线上直线段特征参数，反映了岩石渗透率、参与渗流过程孔隙大小分布、流体粘度、界面张力和润湿性等参数。

这里，将平面径向流相对渗透率模型与一维线性流相对渗透率模型进行了对比，分析了相对渗透率比值与含水饱和度关系曲线上直线段特征参数与饱和度区间，见表 2、图 6。

对比结果表明，水驱过程中，与油藏深处线性渗流不同，油藏近井地带油水两相径向渗流中，线性饱和度区间小。在油藏深处，相对渗透率比值与含水饱和度存在一种“直线关系”，呈线性渗流特征；而在油藏近井地带，这种油水相这种线性渗流持续时间短，油井生产过程中主要以非线性渗流关系为主。

表 2还表明，近井地带渗流特征参数A、B与油藏深处也不相同。油藏深处油水流动可以近似为线性流动，而近井地带为径向流动，越靠近井筒，流动速度越快，油水相互作用，相互干扰越强，因此，渗流特征参数不同。

此外，含水饱和度在60%以上时，在注水井附近，径向流相对渗透比值比线性流低，这说明，该地带水相有着非常强的驱油能力，将更多的剩余油驱赶到油层深处；当含水饱和度低于40%，即在油井周围，径向流相对渗透率比值比线性流高，说明在该地带，油相流动能力更强，更容易在油井中汇集。

分析表明，上述方法是一种行之有效的径向流分析方法，能准确描述两相径向流渗流过程，根据油藏生产动态确定油藏近井地带油水相对渗流关系。

3 结语

依据水驱油实验建立了近井地带油水径向渗流的力学模型，得到了相应的解析解，得到了油藏近井地带油水相对渗透率比值与含水饱和度函数关系，并对近井地带油水径向渗流特点进行了分析。

与Buckley-Leverrt一维渗流模型相比，笔者所建立的径向流渗流模型对油藏近井地带实际发生的径向驱油过程和驱油结果更为贴切。

近井地带油水两相径向流与油藏深处线性渗流特征差异显著，近井地带径向流中，油相渗流与水相渗流存在十分复杂的非线性关系。

径向流渗流模型是一种行之有效的分析方法，为准确描述油水井生产动态，分析油水两相在油藏中渗流特征提供了理论依据。