CN103161436B - 一种稠油热采水平井试井解释方法 - Google Patents

Abstract

本发明提供一种稠油热采水平井试井解释方法，包括：根据稠油水平井的地下流体和热场分布特征，采用时变线源复合试井模型；根据所述时变线源复合试井模型构造线源函数；根据稠油水平井的边界条件选择板源函数；叠加所述线源函数和板源函数，获得稠油热采水平井源函数；由所述稠油热采水平井源函数，获得无井储井底压力函数的拉普拉斯空间解；由所述无井储井底压力函数的拉普拉斯空间解，获得有井储井底压力函数的真实空间解。本发明将热波及区假设为以水平井水平段为轴的圆柱体，将新假设下稠油水平井焖井测试过程和高温生产测试过程的井底压力求解问题归结为一个板源函数和一个线源函数的叠加，为稠油热采水平井试井解释提供了新方法。

Description

一种稠油热采水平井试井解释方法

技术领域

本发明涉及试井技术领域，具体地，涉及一种稠油热采水平井试井解释方法。

背景技术

随着发现大型稀油油藏的可能性不断降底，稠油油藏的高效开发成为我国石油工业未来重要发展方向之一。与稀油油藏相比，稠油油藏开发由于地下原油粘度大、水油流度比高、水驱开发效果差，所以多采用或者辅助采用热力采油、化学采油的开发方式。由于随着温度的增加，稠油粘度在一定范围内下降非常明显，具有很好的粘温相关性，因此，热力采油成为稠油开发过程中必须研究的课题之一。

对于稠油油藏热力采油，水平井开采可以使油藏获得更强的蒸汽注入能力和原油生产能力。同时，水平井在开采薄层稠油油藏时具有直井无法比拟的效果。但目前稠油油藏水平井热采模型常借鉴直井热采模型，简化为平行于水平面的波及区和非波及区构成的复合模型。这种简化只考虑了热在水平面的波及，而没有考虑垂向上热传导形成的复合特征，与实际油藏动用特征差异较大，制约了稠油油藏动态监测成果的应用。

目前稠油油藏高温动态监测工艺已较为成熟，高精度、耐高温压力计技术成熟。为了高效开发稠油油藏，确保各种综合治理方案合理有效，进行水平井焖井和高温生产试井分析方法研究以准确获取稠油油藏动态参数是其中关键，具有重要的现实意义。水平井在稠油油藏热采方面具有直井不可比拟的优势，其中水平井热采试井是目前为数不多的一种广泛应用于稠油油藏的动态监测技术，但稠油油藏水平井热采试井解释面临数学模型建立和求解困难的瓶颈，严重制约稠油油藏水平井热采动态认识能力。

发明内容

本发明实施例的主要目的在于提供一种稠油热采水平井试井解释方法，以解决现有的水平井热采复合模型简化不合理、数学模型求解困难的问题。

为了实现上述目的，本发明实施例提供一种稠油热采水平井试井解释方法，包括：

根据稠油水平井的地下流体和热场分布特征，采用时变线源复合试井模型；

根据所述时变线源复合试井模型构造线源函数；

根据稠油水平井的边界条件选择板源函数；

叠加所述线源函数和板源函数，获得稠油热采水平井源函数；

由所述稠油热采水平井源函数，获得无井储井底压力函数的拉普拉斯空间解；

由所述无井储井底压力函数的拉普拉斯空间解，获得有井储井底压力函数的真实空间解。

借助于上述技术方案，本发明针对水平井热采时地下流体物性的分布特征，将热波及区假设为以水平井水平段为轴的圆柱体，将新假设下稠油水平井焖井测试过程和高温生产测试过程的井底压力求解问题归结为一个板源函数和一个线源函数的叠加，为稠油油藏水平井焖井和高温生产试井科学分析提供了简便、可行的新方法。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例一提供的稠油热采水平井试井解释方法的流程示意图；

图2是本发明实施例一提供的构造线源函数的具体流程示意图；

图3是本发明实施例一提供的获得无井储井底压力函数的真实空间解的具体流程示意图；

图4是本发明实施例一提供的获得无井储井底压力函数的拉普拉斯空间解的具体流程示意图；

图5是本发明实施例一提供的无井储井底压力变化曲线离散化示意图；

图6是本发明实施例一提供的获得有井储井底压力函数的真实空间解的具体流程示意图；

图7是本发明实施例二提供的应用稠油热采水平井试井解释方法对稠油水平井进行试井解释的具体流程图；

图8是本发明实施例三提供的实例井水平井焖井阶段拟合曲线；

图9是本发明实施例四提供的实例井水平井高温生产阶段拟合曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，本发明实施例中所指的“多元热流体”包括稠油热采过程中注入的热水、蒸汽、非凝结气，其中，非凝结气包括N₂和CO₂；本发明实施例中所指的“水平井焖井测试过程”为水平井注入多元热流体后，关井测取关井后井底压力、温度随时间变化关系的试井测试过程；本发明实施例中所指的“水平井高温生产测试过程”为水平井注入多元热流体并焖井后，重新开井以一定产量生产，同时测取井底压力、温度随时间变化关系的试井测试过程；本发明实施例中所指的“热采方式”为蒸汽驱或者蒸汽吞吐。

实施例一

本实施例提供一种稠油热采水平井试井解释方法，如图1所示，该方法包括：

步骤S11，根据稠油水平井的地下流体和热场分布特征，采用时变线源复合试井模型；

步骤S12，根据所述时变线源复合试井模型构造线源函数；

步骤S13，根据稠油水平井的边界条件选择板源函数；

步骤S14，叠加所述线源函数和板源函数，获得稠油热采水平井源函数；

步骤S15，由所述稠油热采水平井源函数，获得无井储井底压力函数的拉普拉斯空间解；

步骤S16，由所述无井储井底压力函数的拉普拉斯空间解，获得有井储井底压力函数的真实空间解。

本实施例步骤S11中，提供以下时变线源复合试井模型，其物理模型假设为：

1)单相微可压缩液体；

2)等温流动；

3)油井半径为r_w，考虑表皮因子S的影响；

4)油井生产前，地层中各点的压力均匀分布，波及区和非波及区分别为P_i1、P_i2；

5)忽略重力和毛管力的影响；

6)线性达西渗流；

7)地层径向复合、等厚、各向同性，井以一产量q生产；

8)地层岩石微可压缩。

其数学计算公式为：

$\left\{\begin{array}{c}\frac{1}{r}\·\frac{\∂}{\∂r}\left(r\frac{\∂{\tilde{P}}_1}{\∂r}\right){|}_{r\≤R}=\frac{\mathrm{\φ}{\mathrm{\μ}}_1{C}_{t1}}{k_1}\frac{\∂{\tilde{P}}_1}{\∂t}\\ \frac{1}{r}\·\frac{\∂}{\∂r}\left(r\frac{\∂{\tilde{P}}_2}{\∂r}\right){|}_{r\≥R}=\frac{\mathrm{\φ}{\mathrm{\μ}}_2{C}_{t2}}{k_2}\frac{\∂{\tilde{P}}_2}{\∂t}\\ {\tilde{P}}_1(r,t){|}_{t=0}={\tilde{P}}_i\\ {\tilde{P}}_2(r,t){|}_{t=0}={\tilde{P}}_i\\ {\tilde{P}}_w={\tilde{P}}_1(r_{\mathrm{we}},t)\\ M_1=2\mathrm{\π}{r}_{\mathrm{we}}h{\left(\frac{\∂{\tilde{P}}_1}{\∂r}\right)}_{r=r_{\mathrm{we}}}\\ \underset{r\→\∞}{\lim }{\tilde{P}}_2(r,t)={\tilde{P}}_i\\ \frac{{\mathrm{\μ}}_1}{{\mathrm{\ρ}}_1{k}_1}\frac{\∂{\tilde{P}}_1}{\∂t}{|}_{r=R}=\frac{{\mathrm{\μ}}_2}{{\mathrm{\ρ}}_2{k}_2}\frac{\∂{\tilde{P}}_2}{\∂t}{|}_{r=R}\\ 2\mathrm{\πrh}\frac{\∂{\tilde{P}}_1}{\∂r}{|}_{r=R}=2\mathrm{\πrh}\frac{\∂{\tilde{P}}_2}{\∂r}{|}_{r=R}\end{array}\right.$ （公式1）

公式1中，各标识符号的具体含义如下：

分别为热/气波及区和非波及区压力对应的拟压力函数；

P₁、P₂分别为热/气波及区和非波及区压力，单位为atm；

ρ₁、ρ₂分别为热/气波及区和非波及区流体密度，单位为g/cm³；

k₁μ₁、k₂μ₂分别为热/气波及区和非波及区流度系数，单位为μm²/mPa.s；

C_t1、C_t2分别为热/气波及区和非波及区综合压缩系数，单位为atm^-1；

为P_i对应的拟压力函数；

R为波及区半径，单位为cm；

M₁为产出质量流，单位为g/s，M₁＝qBρ₁，其中，B为体积系数，无因次。

所述拟压力函数与压力的关系式：

$\tilde{P}={\∫}_0^P\frac{\mathrm{\ρk}}{\mathrm{\μ}}\mathrm{dP}$ （公式2）

公式2中各标示符号的具体含义如下：

为计算压力值对应的拟压力函数；

P为计算压力值，单位为atm；

ρ为探测半径内流体密度，单位为g/cm³；

kμ为探测半径内流度系数，单位为μm²/mPa.s。

优选的，如图2所示，本实施例的步骤S12具体包括：

步骤S121，对所述时变线源复合试井模型进行无因次化和拉普拉斯变换，获得无井储井底拟压力函数的拉普拉斯空间解；

具体的，本步骤中，进行无因次化采用如下公式：

$\left\{\begin{array}{c}{\tilde{P}}_{1D}=\frac{2\mathrm{\πh}({\tilde{P}}_{i1}-{\tilde{P}}_1)}{M_1}\\ {\tilde{P}}_{2D}=\frac{2\mathrm{\πh}({\tilde{P}}_{i2}-{\tilde{P}}_2)}{M_1}\\ t_D=\frac{k_{1}t}{\mathrm{\φ}{\mathrm{\μ}}_1{C}_{t1}{r}_{\mathrm{we}}^2}\\ r_D=\frac{r}{r_{\mathrm{we}}}\end{array}\right.$ （公式3）

进行拉普拉斯变换采用如下公式：

$\stackrel{\‾}{P_D}={\∫}_0^{\∞}{\tilde{P}}_D{e}^{-u{t}_D}d{t}_D$ （公式4）

最终，获得的无井储井底拟压力函数的拉普拉斯空间解如下：

$\left\{\begin{array}{c}\stackrel{\‾}{P_{\mathrm{wD}}}=A_1{I}_0\left(\sqrt{u}\right)+B_1{K}_0\left(\sqrt{u}\right)\\ B_1=\frac{1}{u}\frac{1}{-\frac{{\mathrm{\β}}_1\frac{{\mathrm{\μ}}_2{\mathrm{\ρ}}_1{k}_1}{{\mathrm{\μ}}_1{\mathrm{\ρ}}_2{k}_2}{K}_0\left(R^{\′}{\mathrm{\β}}_2\right)K_1\left(R^{\′}{\mathrm{\β}}_1\right)-{\mathrm{\β}}_2{K}_1\left(R^{\′}{\mathrm{\β}}_2\right)K_0\left(R^{\′}{\mathrm{\β}}_1\right)}{{\mathrm{\β}}_2{K}_1\left(R^{\′}{\mathrm{\β}}_2\right)I_0\left(R^{\′}{\mathrm{\β}}_1\right)+{\mathrm{\β}}_1\frac{{\mathrm{\μ}}_2{\mathrm{\ρ}}_1{k}_1}{{\mathrm{\μ}}_1{\mathrm{\ρ}}_2{k}_2}{K}_0\left(R^{\′}{\mathrm{\β}}_2\right)I_1\left(R^{\′}{\mathrm{\β}}_1\right)}{\mathrm{\β}}_1{I}_1\left({\mathrm{\β}}_1\right)+{\mathrm{\β}}_1{K}_1\left({\mathrm{\β}}_1\right)}\\ A_1=\frac{1}{u}\frac{1}{-{\mathrm{\β}}_1{I}_1\left({\mathrm{\β}}_1\right)+\frac{{\mathrm{\β}}_2{K}_1\left(R^{\′}{\mathrm{\β}}_2\right)I_0\left(R^{\′}{\mathrm{\β}}_1\right)+{\mathrm{\β}}_1\frac{{\mathrm{\μ}}_2{\mathrm{\ρ}}_1{k}_1}{{\mathrm{\μ}}_1{\mathrm{\ρ}}_2{k}_2}{K}_0\left(R^{\′}{\mathrm{\β}}_2\right)I_1\left(R^{\′}{\mathrm{\β}}_1\right)}{{\mathrm{\β}}_1\frac{{\mathrm{\μ}}_2{\mathrm{\ρ}}_1{k}_1}{{\mathrm{\μ}}_1{\mathrm{\ρ}}_2{k}_2}{K}_0\left(R^{\′}{\mathrm{\β}}_2\right)K_1\left(R^{\′}{\mathrm{\β}}_1\right)-{\mathrm{\β}}_2{K}_1\left(R^{\′}{\mathrm{\β}}_2\right)K_0\left(R^{\′}{\mathrm{\β}}_1\right)}{\mathrm{\β}}_1{K}_1\left({\mathrm{\β}}_1\right)}\\ R^{\′}=\frac{R}{r_{\mathrm{we}}}\\ {\mathrm{\β}}_1=\sqrt{u}\\ {\mathrm{\β}}_2=\sqrt{u{k}_1{\mathrm{\μ}}_2{C}_{t2}/\left(k_2{\mathrm{\μ}}_1{C}_{t1}\right)}\end{array}\right.$

(公式5)

公式3、4、5中，各标识符号的具体含义如下：

为波及区无因次拟压力函数；h为油层厚度；M₁为产出质量流；分别为波及区和非波及区压力对应的拟压力函数；分别为波及区压力对应的拟压力函数和非波及区压力对应的拟压力函数；为非波及区无因次拟压力函数；t_D为无因次时间；t为时间；φ为平均孔隙度；k₁为波及区渗透率；μ₁为波及区流体粘度；C_t1为波及区的综合压缩系数；r_D为无因次距离；r为半径；r_we为有效井筒半径；为拉普拉斯变换的像函数；为对时变线源复合试井模型进行无因次化和拉普拉斯变换后，求得的无井储井底拟压力函数的拉普拉斯空间解；I₀为0阶第一类虚宗量贝塞尔函数。

步骤S122，使用Stehfest反演方法，由所述无井储井底拟压力函数的拉普拉斯空间解得到无井储井底拟压力函数的真实空间解；

具体的，该步骤采用如下公式：

$\left\{\begin{array}{c}P(r,t)=\frac{\ln 2}{t}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{V}_i\stackrel{\‾}{P}(r,u)\\ V_i={(-1)}^{\frac{N}{2}+i}\underset{n=\left[\frac{i+1}{2}\right]}{\overset{\min (i,\frac{N}{2})}{\mathrm{\Σ}}}\frac{n^{\frac{N}{2}}\left(2n\right)!}{(\frac{N}{2}-n)!n!(n-1)!(i-n)!(2n-i)!}\end{array}\right.$ （公式6）

公式6中，各标示符号的具体含义为：P(r,t)为无井储井底拟压力函数的真实空间解；为无井储井底拟压力函数的拉普拉斯空间解；N为4－16之间的偶数。

步骤S123，通过将所述无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解，所述时态域是自变量为时间、压力、温度的多维空间，所述时间域是自变量为时间的一维空间；

具体的，如图3所示，步骤S123具体包括：

步骤S1231，对所述无井储井底拟压力函数进行离散化；

该步骤采用如下公式：

$\left\{\begin{array}{c}{\tilde{P}}_n=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\ρ}\left(P_i\right)k\left(P_i\right)}{\mathrm{\μ}\left(P_i\right)}\mathrm{\Δ}{P}_i\\ {\tilde{P}}_n-{\tilde{P}}_{n-1}=\frac{\mathrm{\ρ}\left(P_n\right)k\left(P_n\right)}{\mathrm{\μ}\left(P_n\right)}\mathrm{\Δ}{P}_n\end{array}\right.$ （公式7）

公式7中，各标识符号的具体含义如下：ρ为探测半径内流体密度；k为探测半径内渗透率；μ为探测半径内流体粘度；k/μ为探测半径内流度；P为压力；ΔP为压力差；为拟压力函数；n为时间序号。

步骤S1232，将所述离散化后的无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解；

该步骤采用如下公式：

$P\left(t_n\right)-P\left(t_{n-1}\right)=\frac{\mathrm{\μ}\left(t_n\right)}{\mathrm{\ρ}\left(t_n\right)k\left(t_n\right)}[\tilde{P}\left(t_n\right)-\tilde{P}\left(t_{n_{-1}}\right)]$ （公式8）

公式8中，各标识符号的具体含义如下：ρ为探测半径内流体密度；k为探测半径内渗透率；μ为探测半径内流体粘度；k/μ为探测半径内流度；P为压力；ΔP为压力差；为拟压力函数；n为时间序号；t为时间。

本实施例的步骤S1231和步骤S1232中，探测半径内流体密度ρ、探测半径内流度k/μ由分流量方程和流体饱和度加权平均获得，具体如下：

公式8中包含密度ρ和流度系数kμ的信息，反映了稠油油藏水平井热采过程中干度、蒸汽、非凝结气、多相的影响，本实施例提供了一套不同压力、温度下密度ρ和流度系数kμ的计算方法；

所述流体密度是不同流体密度按饱和度加权平均得到，其中自由变量包括S_o,S_g,S_w,ρ_o,ρ_g,ρ_w,T,P。

其中，S_o为含油饱和度，小数；

S_g为含气饱和度，小数；

S_w为含水饱和度，小数；

ρ_o为油相密度，单位为g/cm³；

ρ_g为气相密度，单位为g/cm³；

ρ_w为水相密度，单位为g/cm³；

T为温度，单位为℃；

P为压力，单位为atm。

含油饱和度S_o按照S_o+S_g+S_w＝1由S_g和S_w计算得到。

含气饱和度S_g采用如下公式，由蒸汽饱和度和非凝结气饱和度求和得到：

S_g＝S_gN2+CO2+S_gsteam（公式9）

公式9中，各标示符号的具体含义为：

S_g为含气饱和度，由气相N₂、CO₂和蒸汽饱和度求和得到，值为小数；

S_gN2+CO2为非凝结气饱和度，值为小数；

S_gsteam为蒸汽饱和度，值为小数；

初始含气饱和度S_g0由初始非凝结气饱和度和初始蒸汽饱和度求和得到。

初始非凝结气饱和度由输入的日注N₂量、日注CO₂量计算求得，初始蒸汽饱和度由输入的日注水量、井底注入蒸汽干度计算求得。

含气饱和度随时间的变化分为非凝结气饱和度随时间的变化和蒸汽饱和度随时间的变化。

所述非凝结气饱和度随时间的变化指：S_gN2+CO2＝S_giN2+CO2。

其中，S_giN2+CO2为多元热流体注入结束时N₂和CO₂的饱和度，值为分数。

所述蒸汽饱和度随时间的变化指：

其中，S_gisteam为多元热流体注入结束时蒸汽的饱和度，值为分数；Δt为测试时间，单位为s；t为蒸汽腔持续时间，单位为s。

优选的，本实施例中，在步骤S123中，若确定试井测试过程为高温生产测试过程，则在不考虑蒸汽腔影响。

具体的，对于水平井高温生产测试过程，不需要考虑蒸汽腔影响，即蒸汽饱和度S_gsteam=0。

初始含水饱和度S_w0由输入的日注水量、井底注入蒸汽干度计算求得。

含水饱和度S_w由分流量方程按照常规方法计算得到。

所述分流量方程由改进的科瓦尔方法建立，具体计算公式如下：

$\left\{\begin{array}{c}{f}_s=\frac{K(1+G)-{[(1+G)K/V_{\mathrm{pi}}]}^{1/2}(1-S_{\mathrm{ot}})}{K-1}\\ K=H\×F\×E=H\×F\×{[0.78+0.22{\left(\frac{{\mathrm{\μ}}_o}{{\mathrm{\μ}}_s}\right)}^{\frac{1}{4}}]}^4\\ F=0.565{\log }_{10}\left(\frac{t_h}{t_v}\right)+0.870\\ \frac{t_h}{t_v}=C\left(\frac{k_v}{k_h}\right)\left(\frac{A}{h}\right)\left(k_{h}h\right)\frac{\mathrm{\Δ\ρ}}{\mathrm{Q\μ}}\\ G=\frac{k_h{k}_{\mathrm{ro}}\mathrm{Ag}}{{\mathrm{\μ}}_o{Q}_{t}B}({\mathrm{\ρ}}_s-{\mathrm{\ρ}}_o)\frac{\sin \mathrm{\β}}{\cos (\mathrm{\α}-\mathrm{\β})}\end{array}\right.$ （公式10）

公式10中，各标识符号的含义如下：

μ_o为地层流体平均粘度，单位为mPa.s；

K_ro为地层流体相渗，值为分数；

g为重力加速度，单位为cm/s²；

f_s为含蒸汽、非凝结气、水率，小数；

表示水平与垂向流体交换能力的比值；

G为重力项；

V_pi为注入蒸汽、非凝结气、热水后的实际孔隙体积，单位为cm³；

S_ot为含油饱和度，小数；

H为非均质系数，log₁₀H＝[V_DP/(1-V_DP)^0.2]，均质系统等于1；

V_DP为戴克斯特拉-帕森斯系数，基于K-h（岩心块数）统计结果得到，并且 $V_{\mathrm{DP}}=\frac{V_{50\%}-V_{84.1\%}}{V_{50\%}}$

E为科瓦尔有效粘度比， $E={[0.78+0.22{\left(\frac{{\mathrm{\μ}}_o}{{\mathrm{\μ}}_s}\right)}^{\frac{1}{4}}]}^4;$

μ_s为饱和原油粘度，单位为mPa.s；

C为井网常数，五点法为2.5271，或线性驱为2.1257；

k_v为垂直渗透率，单位为μm²；

k_h为横向渗透率，单位为μm²；

A为井网规模，英亩；

Δρ为蒸汽、非凝结气、水混合流体与原油密度差，单位为g/cm³；

Q为注入速度，单位为桶/天；

μ为油、气、水混合粘度，单位为mPa.s，由常规方法加权得到；

α为地层倾角，弧度；

β为界面倾角，弧度；

Q_t为总生产率；

k_H为横向渗透率，单位为μm²；

ρ_s为蒸汽、非凝结气、水混合密度，单位为g/cm³，由常规方法加权得到。

由分流量方程f_s可计算得到蒸汽、非凝结气、水的总饱和度S_s。

采用如下公式，根据注入蒸汽、非凝结气与注入水量的比例，可计算得到含水饱和度S_w：

$S_w=\frac{q_w}{q_w+q_g}{S}_s$ （公式11）

公式11中，q_w为注入水的速度，单位为cm³/s；q_g为注入气体（蒸汽、非凝结气）的速度，单位为cm³/s。

油相密度ρ_o、气相密度ρ_g、水相密度ρ_w与温度T和压力P相关，具体关系式由实验室测得或使用已有经验关系式。其中，根据经验构建温度T的热计算式为：

T＝Tt^a+1（公式12）

公式12中，T为测试开始时温度，单位为℃；a为系数，通过拟合测试温度得到。

根据经验选取原始地层压力作为典型压力：P_i+1＝P_i

所述流度系数计算公式如下：

$\left\{\begin{array}{c}\frac{k}{\mathrm{\μ}}=k_e(\frac{k_{\mathrm{ro}}}{{\mathrm{\μ}}_o}+\frac{k_{\mathrm{rw}}}{{\mathrm{\μ}}_w}+\frac{k_{\mathrm{rg}}}{{\mathrm{\μ}}_g})\\ k_{\mathrm{rw}}=k_{\mathrm{rw}}\left(S_w\right)\\ k_{\mathrm{rg}}=k_{\mathrm{rg}}\left(S_g\right)\\ k_{\mathrm{ro}}=k_{\mathrm{ro}}\left(S_o\right)\end{array}\right.$ （公式13）

公式13中，各标识符号的具体含义如下：

k_e为绝对渗透率，单位为μm²；

k_rw为水相相对渗透率，单位为μm²；

k_ro为油相相对渗透率，单位为μm²；

k_rg为气相相对渗透率，单位为μm²；

μ_o为油相粘度，单位为mPa.s；

μ_w为水相粘度，单位为mPa.s；

μ_g为气相粘度，单位为mPa.s。

对于地层原油粘度，当地层压力P不大于饱和压力P_b时，使用Beggs和Robinson（1975）公式描述：

${\mathrm{\μ}}_o=\frac{10.715}{{(5.615R_s+100)}^{0.515}}{{\mathrm{\μ}}_{\mathrm{oD}}}^{\frac{5.44}{{(5.615R_s+150)}^{0.338}}}$ （公式14）

公式14中，μ_oD为脱气原油粘度，单位为mPa.s；R_s为溶解气油比，单位为cm³/cm³。

当地层压力P大于泡点压力P_b时，此时不饱和原油粘度表示为：

μ_o(P)＝μ_o(P_b)+A'(P-P_b)（公式15）

公式15中，A'为P＞P_b时原油粘度与压力曲线斜率。

参考Marhoun方法（1988），溶解度与压力和温度的关系式：

$R_s=A_1{P}^{B_1}{T}^{C_1}$ （公式16）

公式16中A₁、B₁、C₁分别为溶解度与压力和温度关系式系数。以某种油样为例，天然气系数分别为：{102819，1.398441，-1.85513}，CO2系数分别为：{312819，1.398441，-1.85513}，N2系数分别为：{0.343518，1.398441，0.2}。

溶解CO₂、N₂对原油粘度的影响参照天然气，使用Beggs和Robinson（1975）公式描述：

$\left\{\begin{array}{c}{\mathrm{\μ}}_o=A_2{{\mathrm{\μ}}_{\mathrm{uD}}}^{B_2}\\ A_2=\frac{b^c}{{(a{R}_s+b)}^c}\\ B_2=\frac{f^g}{{(e{R}_s+f)}^g}\end{array}\right.$ （公式17）

公式17中，a、b、c、e、f、g分别为粘度关系式系数。以某种油样为例，天然气系数分别为：{5.615，100，0.515，5.615，150，0.338}，CO₂系数分别为：{5.615，100，0.515，5.615，150，0.338}，N₂系数分别为：{1.615，100，0.515，1.615，150，0.338}。

步骤S124，对所述无井储井底压力函数的真实空间解求导，获得线源函数。

该步骤中，无井储井底压力函数的真实空间解为：

$P_i-P_w=\mathrm{\ΔP}(M,t)=\frac{q}{\mathrm{\φcL}}{\∫}_0^t{G}_{\mathrm{\λ}}(M,\mathrm{\τ})\mathrm{d\τ}$ （公式18）

对公式18的两端进行求导，转换为径向坐标系：

$\frac{\∂\mathrm{\ΔP}(r,t)}{\∂t}=\frac{q}{\mathrm{\φcL}}{G}_{\mathrm{\λ}}(r,t)$ （公式19）

得到的线源函数为：

$G_{\mathrm{\λ}}(r,t)=\frac{\mathrm{\φcL}}{q}\frac{\∂\mathrm{\ΔP}(r,t)}{\∂t}$ （公式20）

公式18、19、20中各标识符号的具体含义如下：

P_i为原始地层压力，单位为atm；

P_w为M点在t时刻无井储井底压力函数的真实空间解，单位为atm；

ΔP(M,t)为M点在t时刻的压力变化，单位为atm；

M为空间中一点的坐标，M＝M(r,X)；

q为产量，单位为cm³/s；

φ为平均孔隙度，单位为小数；

c为综合压缩系数，单位为atm^-1；

L为水平井水平段半长，单位为cm；

t为时间，单位为s；

G_λ为所求的线源函数；

τ为被积时间变量，单位为s；

X为一维方向；

r为M点在极坐标极轴上的长度，单位为cm。

本实施例中，步骤S13具体用于首先根据地质条件判断稠油水平井的边界条件（例如为封闭边界、定压边界、混合边界），然后对照现有文献提供的板源函数表选择满足稠油水平井边界条件的板源函数，现有文献提供了封闭边界、定压边界、混合边界分别对应的无限大板源函数，因此本步骤可根据现有文献提供的板源函数中选择。

优选的，本实施例步骤S14中，叠加所述线源函数和板源函数，具体为：利用Newman乘积方法将所述线源函数和板源函数叠加。

具体的，该步骤采用如下公式：

S(M,t)＝G_λ(r,t)×G_X(X,t)（公式21）

公式21中，G_X(X,t)为板源函数；G_λ(r,t)为线源函数，S(M,t)为M处的稠油热采水平井源函数。

优选的，如图4所示，本实施例中步骤S15具体包括：

步骤S151，应用格林函数法对所述稠油热采水平井源函数进行积分，获得无井储井底压力函数的真实空间解；

具体的，该步骤应用格林函数法，对叠加得到的稠油热采水平井源函数S(M,t)进行积分，即可获得无井储井底压力变化曲线ΔP(M,t)：

$\mathrm{\ΔP}(M,t)=P_i-P_w=\frac{q}{\mathrm{\φcL}}{\∫}_0^{t}S(M,\mathrm{\τ})\mathrm{d\τ}$ （公式22）

公式22中，各标识符号的具体含义如下：

P_w为水平井开采下M点在t时刻无井储井底压力函数的真实空间解，单位为atm；

P_i为原始地层压力，单位为atm；

ΔP(M,t)为水平井开采下M点在t时刻的压力变化，单位为atm；

q为产量，单位为cm³/s；

φ为平均孔隙度，单位为小数；

c为综合压缩系数，单位为atm^-1；

L为水平井水平段半长，单位为cm；

t为时间，单位为s；

τ为被积时间变量，单位为s。

步骤S152，对所述无井储井底压力函数的真实空间解进行离散化处理，得到离散的无井储井底压力函数的真实空间解；

具体的，如图5所示，该步骤中将无井储井底压力变化曲线离散化成若干段，各段压力变化分别为：dP₁,dP₂,dP₃…,dP_n

则无井储井底压力函数的真实空间解可以表示为：

$\left\{\begin{array}{c}{P}_w=\mathrm{\&theta;}(t-t_0)d{P}_1+\mathrm{\&theta;}(t-t_1)d{P}_2+\mathrm{\&theta;}(t-t_2)d{P}_3+\&CenterDot;\&CenterDot;\&CenterDot;+\mathrm{\&theta;}(t-t_{n-1})d{P}_n\\ \mathrm{\&theta;}(t-t_i)=\left\{\begin{array}{cc}1& t\&GreaterEqual;t_i\\ 0& t<t_i\end{array}\right.\end{array}\right.$ （公式23）

公式23中，各标识符号的具体含义如下：

P_w为水平井开采下M点在t时刻的无井储井底压力函数的真实空间解，单位为atm；

θ(t-t_i)为阶跃函数，i＝1,2,…,n；

dP_i为t_i时刻对应P_w减去t_i-1时刻对应P_w的差，单位为atm；

t为时间，单位为s；

t_i为时间离散点，单位为s；

步骤S153，对所述离散的无井储井底压力函数的真实空间解进行拉普拉斯数值变换，获得无井储井底压力函数的拉普拉斯空间解。

具体的，该步骤中使用如下公式：

$\stackrel{\&OverBar;}{P}=e^{-u{t}_0}\frac{d{P}_1}{u}+e^{-u{t}_1}\frac{d{P}_2}{u}+e^{-u{t}_2}\frac{d{P}_3}{u}+\&CenterDot;\&CenterDot;\&CenterDot;+e^{-u{t}_{n-1}}\frac{d{P}_n}{u}$ （公式24）

公式24中，各标识符号的具体含义如下：

为无井储井底压力函数的拉普拉斯空间解；

e为自然常数；

u为拉普拉斯变量；

t_i为时间离散点，单位为s；

dP_i为t_i时刻对应P_w减去t_i-1时刻对应P_w的差值，单位为atm，i＝1,2,…,n。

优选的，如图6所示，本实施例中步骤S16具体包括：

步骤S161，应用杜哈美原理，将所述无井储井底压力函数的拉普拉斯空间解转换为有井储井底压力函数的拉普拉斯空间解；

具体的，该步骤采用如下公式：

$\stackrel{\&OverBar;}{P_{\mathrm{wD}}\left(C_D\right)}=\frac{\stackrel{\&OverBar;}{P_{\mathrm{wD}}(C_D=0)}}{1+C_D{u}^2\stackrel{\&OverBar;}{P_{\mathrm{wD}}(C_D=0)}}$ （公式25）

公式25中，各标示符号的具体含义如下：

为有井储井底拟压力函数的拉普拉斯空间解；

C_D为无因次化后井储系数，无因次定义式为：C为井储系数，单位为cm³/atm。

步骤S162，对所述有井储井底压力函数的拉普拉斯空间解进行Stehfest数值反演，获得有井储井底压力函数的真实空间解；

具体的，该步骤采用如下公式：

$\left\{\begin{array}{c}P(r,t)=\frac{\ln 2}{t}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_i\stackrel{\&OverBar;}{P}(r,u)\\ V_i={(-1)}^{\frac{N}{2}+i}\underset{n=\left[\frac{i+1}{2}\right]}{\overset{\min (i,\frac{N}{2})}{\mathrm{\&Sigma;}}}\frac{n^{\frac{N}{2}}\left(2n\right)!}{(\frac{N}{2}-n)!n!(n-1)!(i-n)!(2n-i)!}\end{array}\right.$ （公式26）

公式26中，各标示符号的具体含义为：P(r,t)为有井储井底压力函数的真实空间解；为有井储井底压力函数的拉普拉斯空间解；N为4－16之间的偶数。

优选的，本实施例中，在步骤S16获得有井储井底压力函数的真实空间解时，若确定试井测试过程为焖井测试过程，则还需采用如下公式对获得的有井储井底压力函数的真实空间解进行叠加：

P_i-P_w＝ΔP_w(-q,t_p+Δt)+ΔP_w(q,Δt)（公式27）

公式27中，P_i为原始地层压力；P_w为焖井测试过程中的井底压力；ΔP_w为焖井测试过程造成的井底压力降落；t_p为焖井测试过程前的生产时间；Δt为焖井测试过程的时间；q为产量。

实施例二

本实施例提供一按照实施例一中稠油热采水平井试井解释方法对稠油水平井进行试井解释的具体应用过程，如图7所示，包括：

步骤S71，输入基础数据；

步骤S72，根据实施例一中的稠油热采水平井试井解释方法，通过调整关键参数，计算得到理论压力和理论压力导数曲线；

步骤S73，利用理论压力和理论压力导数曲线，拟合实测压力和压力导数曲线；

步骤S74，根据拟合结果得到油藏参数。

其中，基础数据包括测试井地质分层、测井解释数据、多元热流体参数、流体物性、测试工艺、原始地层压力、原油泡点压力、管柱结构、试井测试数据等。

所述地质分层数据包括测试井前期地质研究成果中的分层数据，用于判断测试目标层的深度。

所述测井解释数据一般包括测井解释的有效厚度、孔隙度、渗透率、非有效储层的厚度和岩性，作为参考和对比的基础；

所述多元热流体参数包括溶解气溶解参数、N₂溶解参数、CO₂溶解参数、溶解气降粘系数、N₂降粘系数、CO₂降粘系数。

所述流体物性包括含气原油粘温曲线、束缚水饱和度、水相最大相对渗透率、残余油饱和度、气体最大相渗、油气水三相相对渗透率曲线。

所述测试工艺包括前期折算注采时间、日注N₂量、日注CO₂量、日注水量（或当量）、井底注入蒸汽干度、井底温度、焖井时间、蒸汽腔持续时间、焖井期间井底温度变化系数、稳定日产油量、稳定日产水量、生产期间井底温度变化系数、产气延续时间、产气量变化系数、CO₂平均浓度、N2平均浓度、测试过程（焖井或高温生产）。

所述管柱结构包括油管直径、套管直径、封隔器位置、油管下入深度、测试仪器深度、测试期间管柱调整顺序。

所述试井测试数据包括试井实测数据（时间、压力、温度数据）、中值滤波后的试井测试数据。

所述关键参数包括井储系数、内区渗透率、表皮因子、加热半径、外区渗透率、井底蒸汽干度、拟稳态压力上升速度。

实施例三

本实施例为按照实施例二提供的具体应用过程对一水平井焖井测试过程进行试井解释的应用实例。本实施例中，选择时变线源复合模型，实例井水平井焖井阶段拟合曲线如图8所示，得到解释数据如表1所示：

表1

内区渗透率，×10-3μm2	2000
		外区渗透率，×10-3μm2	2000
吸气段长度，m	40
		边界距离，m	无响应
井底蒸汽干度，f	0.003
		加热半径，m	2.5
井储系数，m3/MPa	10
		表皮因子	0.1

实施例四

本实施例为按照实施例二提供的具体应用过程对一水平井高温生产测试过程进行试井解释的应用实例。本实施例中，选择时变线源复合模型，实例井水平井高温生产阶段拟合曲线如图9所示，得到解释数据如表2所示：

表2

内区渗透率，×10-3μm2	2000
		外区渗透率，×10-3μm2	2000
产液段长度，m	40
		边界距离，m	无响应
加热半径，m	1
		井储系数，m3/MPa	1
表皮因子	0.1

综上所述，本发明实施例提供的稠油热采水平井试井解释方法具有以下有益效果：

（1）本发明针对水平井热采时地下流体物性的分布特征，提出了一种更合理的物理模型，与现有的水平井试井解释方法主要参考直井假设的情况不同，将热波及区假设为以水平井水平段为轴的圆柱体，更符合水平井热采时热的实际地下分布情况；

（2）本发明提出了一种新的稠油热采水平井源函数构建方法，通过选择符合边界条件的板源函数和构建具有分区特征的复合模型线源函数，利用Newman乘积方法，实现了水平井热采井底压力解的计算；构建不同边界条件下的源函数，可快速、方便求解不同边界条件下的水平井热采井底压力；

（3）本发明针对水平井热采条件下干度、蒸汽、非凝结气、多相的影响，引入拟压力函数时态转化方法，采用拉普拉斯变换、Stehfest数值反演、杜哈美原理等方法，建立了分别针对水平井焖井和水平井高温生产井底压力的计算方法，提出了一套对应的试井解释新方法。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种稠油热采水平井试井解释方法，其特征在于，包括：

根据稠油水平井的地下流体和热场分布特征，采用时变线源复合试井模型；

根据所述时变线源复合试井模型构造线源函数；

根据稠油水平井的边界条件选择板源函数；

叠加所述线源函数和板源函数，获得稠油热采水平井源函数；

由所述稠油热采水平井源函数，获得无井储井底压力函数的拉普拉斯空间解；

由所述无井储井底压力函数的拉普拉斯空间解，获得有井储井底压力函数的真实空间解；

其中，所述时变线源复合试井模型为：

1)单相微可压缩液体；

2)等温流动；

3)油井半径为r_w，考虑表皮因子S的影响；

4)油井生产前，地层中各点的压力均匀分布，波及区和非波及区分别为P_i1、P_i2；

5)忽略重力和毛管力的影响；

6)线性达西渗流；

7)地层径向复合、等厚、各向同性，井以一产量q生产；

8)地层岩石微可压缩；

其数学计算公式为：

$\left\{\begin{array}{c}\frac{1}{r}\&CenterDot;\frac{\&PartialD;}{\&PartialD;r}(r-\frac{\&PartialD;{\tilde{P}}_1}{\&PartialD;r}){|}_{r\&le;R}=\frac{\mathrm{\&phi;}{\mathrm{\&mu;}}_1{C}_{t1}}{k_1}\frac{\&PartialD;{\tilde{P}}_1}{\&PartialD;t}\\ \frac{1}{r}\&CenterDot;\frac{\&PartialD;}{\&PartialD;r}\left(r\frac{\&PartialD;{\tilde{P}}_2}{\&PartialD;r}\right){|}_{r\&GreaterEqual;R}=\frac{\mathrm{\&phi;}{\mathrm{\&mu;}}_2{C}_{t2}}{k_2}\frac{\&PartialD;{\tilde{P}}_2}{\&PartialD;t}\\ {\tilde{P}}_1(r,t){|}_{t=0}={\tilde{P}}_i\\ {\tilde{P}}_2(r,t){|}_{t=0}={\tilde{P}}_i\\ {\tilde{P}}_w={\tilde{P}}_1(r_{\mathrm{we}},t)\\ M_1=2\mathrm{\&pi;}{r}_{\mathrm{we}}h{\left(\frac{\&PartialD;{\tilde{P}}_1}{\&PartialD;r}\right)}_{r=r_{\mathrm{we}}}\\ \underset{r\&RightArrow;\&infin;}{\mathrm{lin}}{\tilde{P}}_2(r,t)={\tilde{P}}_i\\ \frac{{\mathrm{\&mu;}}_1}{{\mathrm{\&rho;}}_1{k}_1}\frac{\&PartialD;{\tilde{P}}_1}{\&PartialD;t}{|}_{r=R}=\frac{{\mathrm{\&mu;}}_2}{{\mathrm{\&rho;}}_2{k}_2}\frac{\&PartialD;{\tilde{P}}_2}{\&PartialD;t}{|}_{r=R}\\ 2\mathrm{\&pi;rh}\frac{\&PartialD;{\tilde{P}}_1}{\&PartialD;r}{|}_{r=R}=2\mathrm{\&pi;rh}\frac{\&PartialD;{\tilde{P}}_2}{\&PartialD;r}{|}_{r=R}\end{array}\right.$

以上公式中，各标识符号的具体含义如下：

分别为热/气波及区和非波及区压力对应的拟压力函数；

P₁、P₂分别为热/气波及区和非波及区压力，单位为atm；

ρ₁、ρ₂分别为热/气波及区和非波及区流体密度，单位为g/cm³；

k₁/μ₁、k₂/μ₂分别为热/气波及区和非波及区流度系数，单位为μm²/mPa.s；

C_t1、C_t2分别为热/气波及区和非波及区综合压缩系数，单位为atm^-1；

为P_i对应的拟压力函数；

R为波及区半径，单位为cm；

M₁为产出质量流，单位为g/s，M₁＝qBρ₁，其中，B为体积系数，无因次；

所述拟压力函数与压力的关系式：

$\tilde{P}={\&Integral;}_0^P\frac{\mathrm{\&rho;k}}{\mathrm{\&mu;}}\mathrm{dP}$

以上公式中各标示符号的具体含义如下：

为计算压力值对应的拟压力函数；

P为计算压力值，单位为atm；

ρ为探测半径内流体密度，单位为g/cm³；

k/μ为探测半径内流度系数，单位为μm²/mPa.s；

其中，根据所述时变线源复合试井模型构造线源函数，包括：

对所述时变线源复合试井模型进行无因次化和拉普拉斯变换，获得无井储井底拟压力函数的拉普拉斯空间解；

使用Stehfest反演方法，由所述无井储井底拟压力函数的拉普拉斯空间解得到无井储井底拟压力函数的真实空间解；

通过将所述无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解，所述时态域是自变量为时间、压力、温度的多维空间，所述时间域是自变量为时间的一维空间；

对所述无井储井底压力函数的真实空间解求导，获得线源函数。

2.根据权利要求1所述的稠油热采水平井试井解释方法，其特征在于，通过将所述无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解，具体包括：

对所述无井储井底拟压力函数进行离散化；

将所述离散化后的无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解。

3.根据权利要求2所述的稠油热采水平井试井解释方法，其特征在于，对所述无井储井底拟压力函数进行离散化时，采用如下公式：

$\left\{\begin{array}{c}{\tilde{P}}_n=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&rho;}\left(P_i\right)k\left(P_i\right)}{\mathrm{\&mu;}\left(P_i\right)}{\mathrm{\&Delta;P}}_i\\ {\tilde{P}}_n-{\tilde{P}}_{n-1}=\frac{\mathrm{\&rho;}\left(P_n\right)k\left(P_n\right)}{\mathrm{\&mu;}\left(P_n\right)}{\mathrm{\&Delta;P}}_n\end{array}\right.$

其中，ρ为探测半径内流体密度；k为探测半径内渗透率；μ为探测半径内流体粘度；k/μ为探测半径内流度；P为压力；ΔP为压力差；为拟压力函数；n为时间序号。

4.根据权利要求3所述的稠油热采水平井试井解释方法，其特征在于，将所述离散化后的无井储井底拟压力函数由时态域转换至时间域时，采用如下公式：

$P\left(t_n\right)-P\left(t_{n-1}\right)=\frac{\mathrm{\&mu;}\left(t_n\right)}{\mathrm{\&rho;}\left(t_n\right)k\left(t_n\right)}[\tilde{P}\left(t_n\right)-\tilde{P}\left(t_{n-1}\right)]$

其中，ρ为探测半径内流体密度；k为探测半径内渗透率；μ为探测半径内流体粘度；k/μ为探测半径内流度；P为压力；ΔP为压力差；为拟压力函数；n为时间序号；t为时间。

5.根据权利要求3或4所述的稠油热采水平井试井解释方法，其特征在于，所述探测半径内流体密度、探测半径内流度由分流量方程和流体饱和度加权平均获得。

6.根据权利要求1所述的稠油热采水平井试井解释方法，其特征在于，叠加所述线源函数和板源函数，具体为：利用Newman乘积方法将所述线源函数和板源函数叠加。

7.根据权利要求1所述的稠油热采水平井试井解释方法，其特征在于，由所述稠油热采水平井源函数，获得无井储井底压力函数的拉普拉斯空间解，具体包括：

应用格林函数法对所述稠油热采水平井源函数进行积分，获得无井储井底压力函数的真实空间解；

对所述无井储井底压力函数的真实空间解进行离散化处理，得到离散的无井储井底压力函数的真实空间解；

对所述离散的无井储井底压力函数的真实空间解进行拉普拉斯数值变换，获得无井储井底压力函数的拉普拉斯空间解。

8.根据权利要求7所述的稠油热采水平井试井解释方法，其特征在于，由所述无井储井底压力函数的拉普拉斯空间解，获得有井储井底压力函数的真实空间解，具体包括：

应用杜哈美原理，将所述无井储井底压力函数的拉普拉斯空间解转换为有井储井底压力函数的拉普拉斯空间解；

对所述有井储井底压力函数的拉普拉斯空间解进行Stehfest数值反演，获得有井储井底压力函数的真实空间解。

9.根据权利要求1所述的稠油热采水平井试井解释方法，其特征在于，获得有井储井底压力函数的真实空间解时，还包括：

若确定试井测试过程为焖井测试过程，则采用如下公式对所述有井储井底压力函数的真实空间解进行叠加：

P_i-P_w＝ΔP_w(-q,t_p+Δt)+ΔP_w(q,Δt)

其中，P_i为原始地层压力；P_w为焖井测试过程中的井底压力；ΔP_w为焖井测试过程造成的井底压力降落；t_p为焖井测试过程前的生产时间；Δt为焖井测试过程的时间；q为产量。

10.根据权利要求1所述的稠油热采水平井试井解释方法，其特征在于，通过将所述无井储井底拟压力函数由时态域转换至时间域，使所述无井储井底拟压力函数的真实空间解转换为无井储井底压力函数的真实空间解时，还包括：

若确定试井测试过程为高温生产测试过程，则不考虑蒸汽腔影响。