CN105089639A - 一种煤层气井井筒流动动态预测方法 - Google Patents

Abstract

本发明公开了一种煤层气井井筒流动动态预测方法，包括：对井筒油管中的液体采样以获取井筒中的真实固体含量；判别井筒中气液固三相流的流动形态，确定气液固三相流的建模参数；依据所述真实固体含量以及所述建模参数，结合动液面位置和井身结构数据，建立井筒内的气液固三相流压力模型；依据所述气液固三相流压力模型，所述真实固体含量以及所述建模参数，结合井身传热参数，建立井筒内的温度分布模型；根据所述压力模型和温度分布模型获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。本发明的方法可以根据井口产液量和产气量，以及套压等物性参数预测任意井型的井筒内任意位置的流态，流速，压力，温度分布等物性参数。

Description

一种煤层气井井筒流动动态预测方法

技术领域

本发明涉及石油工程技术领域，具体地说，涉及一种煤层气井井筒流动动态预测方法。

背景技术

目前，煤层气排采过程中多采用油管产水套管产气的生产方式，与常规天然气生产方式差别较大。同时由于煤质具有较脆、胶结性差、易碎和易坍塌等特点，前期压裂及排采过程中的生产压差和流体作用会造成煤层破坏产生煤粉。由于产出的煤粉中部分颗粒极小，同时煤基质呈现多孔结构，本身视密度较小，因此部分煤粉能够紧随井筒内的流体流动，这也使得煤层气井井筒套管中的流动呈现为煤层气、水、煤粉构成的气液固三相流。

现有关于井筒内三相流的研究主要针对油气水三相流。而对于气液固构成的三相流，由于在油气开采过程中较少遇到，同时固相的存在使得已有的较为完备的气液两相流流型判别方法不再适用，因此目前关于井筒内气液固三相流流动特征的研究极为少见。

基于上述情况，亟需一种煤层气井井筒流动态预测方法准确地预测和分析沿井筒流体压力、温度、相态等的变化特征。

发明内容

本发明针对现有技术的不足，提出了一种煤层气井井筒流动动态预测方法，包括：

对井筒油管中的液体采样以获取井筒中的真实固体含量H_s；测量井口的液相体积流量和气相体积流量，以及井口温度和井口压力；

判别井筒中气液固三相流的流动形态，并根据井筒中的真实固体含量H_s确定气液固三相流的建模参数；

依据所述真实固体含量H_s以及所述建模参数，结合动液面位置和井身结构数据，基于贝格斯-布里尔方法建立井筒内的气液固三相流压力模型；

依据所述气液固三相流压力模型，所述真实固体含量H_s以及所述建模参数，结合井身传热参数，根据能量守恒性质建立井筒内的温度分布模型；

从井口开始对井筒依次划分为若干连续的子井段，将测量得到的井口的液相体积流量和气相体积流量，以及井口温度和井口压力做为初始值，在所述连续的子井段内根据所述压力模型和温度分布模型耦合迭代计算直至井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

根据本发明的一个实施例，还包括根据所述压力和温度沿井筒深度的分布结果以及气体状态方程获得井筒内气相密度沿井筒深度的分布结果。

根据本发明的一个实施例，所述判别井筒中气液固三相流的流动形态包括

由井筒中的液相体积流量参数和气相体积流量参数确定气液固混合物的平均流速 $v_m=\frac{Q_l+Q_g}{0.25\mathrm{\π}{D}^2}$ 和无滑脱持液率 $E_L=\frac{Q_l}{Q_l+Q_g};$

根据弗洛德准数和无滑脱持液率应用贝格斯-布里尔方法划分井筒中气液固三相流的流动形态；

其中，Q_l为井筒中的液相体积流量参数，Q_g为井筒中的气相体积流量参数，g为重力加速度，D为管道直径。

根据本发明的一个实施例，所述根据真实固体含量H_s确定气液固三相流的建模参数包括根据真实固体含量H_s和井筒中气液固三相流的流动形态的流型参数计算真实液体含量H_l(θ)＝(1-H_s)H'_l(θ)，真实气体含量H_g(θ)＝(1-H_s)[1-H'_l(θ)]，

其中，H'_l(θ)＝H_l(0)ψ，H'_l(θ)为倾角为θ的气液两相流动的液体含量，H_l(0)为同样流型参数下水平流的液体含量，ψ为倾斜校正系数，θ为井筒管道与水平方向的夹角。

根据本发明的一个实施例，所述根据真实固体含量H_s确定气液固三相流的建模参数还包括

根据真实固体含量H_s，真实液体含量H_l(θ)和无滑脱持液率E_L计算气液固三相流沿程阻力系数λ＝λ′·e^s；其中，

无滑脱的沿程阻力系数 ${\mathrm{\λ}}^{\′}={\left[21g\left(\frac{R_e^{\′}}{4.52231g{R}_e^{\′}-3.8215}\right)\right]}^{-2};$

R_e′为无滑脱的雷诺数，由含有固相颗粒的液相粘度μ_ls＝μ_l(1+2.5H_s)确定，μ_l为液相粘度；

指数 $s=\frac{\ln Y}{-0.0523+3.182\ln Y-0.8725{\left(\ln Y\right)}^2+0.01853{\left(\ln Y\right)}^4},$

其中，E_L为无滑脱持液率，H_l(θ)为真实液体含量，θ为井筒管道与水平方向的夹角。

根据本发明的一个实施例，所述气液固三相流压力模型由下式表示：

$\frac{\mathrm{dp}}{\mathrm{dz}}-\frac{[{\mathrm{\ρ}}_l{H}_l\left(\mathrm{\θ}\right)+{\mathrm{\ρ}}_g{H}_g\left(\mathrm{\θ}\right)+{\mathrm{\ρ}}_s{H}_s]g\sin \mathrm{\θ}+\frac{\mathrm{\λ}{\mathrm{Gv}}_m}{2\mathrm{DA}}}{1-\frac{[{\mathrm{\ρ}}_l{H}_l\left(\mathrm{\θ}\right)+{\mathrm{\ρ}}_g{H}_g\left(\mathrm{\θ}\right)+{\mathrm{\ρ}}_s{H}_s]v_m{v}_{\mathrm{sg}}}{p}},$

其中，ρ_l为液相密度，ρ_g为气相密度，ρ_s为固相密度，p为气液固三相混合物的压力，z为沿井筒轴向流动的距离，g为重力加速度，G为气液固三相混合物的质量流量，A为管道横截面积，D为管道直径，v_m为气液固三相混合物的平均流速，v_sg为气相表观流速；H_s为真实固体含量，H_l(θ)为真实液体含量，H_g(θ)为真实气体含量，θ为井筒管道与水平方向的夹角；λ为沿程阻力系数。

根据本发明的一个实施例，所述温度分布模型由下式表示：

$C_{\mathrm{pm}}\frac{\mathrm{dT}}{\mathrm{dz}}-C_{\mathrm{pm}}{C}_{\mathrm{Jm}}\frac{\mathrm{dp}}{\mathrm{dz}}+\frac{v_{m}d{v}_m}{\mathrm{dz}}+g\sin \mathrm{\θ}+\frac{\mathrm{\λ}{v}_m^2}{2D}=-\frac{\mathrm{dq}}{\mathrm{dz}}$

其中，由所述气液固三相流压力模型确定，T为井筒温度，C_pm为气液固三相混合物的平均定压比热容，C_Jm为气液固三相混合物的焦耳—汤姆逊数，q为径向热流量，p为气液固三相混合物的压力，v_m为气液固三相混合物的平均流速，λ为气液固三相流沿程阻力系数，θ为井筒管道与水平方向的夹角，D为管道直径，g为重力加速度，z为沿井筒轴向流动的距离。

根据本发明的一个实施例，所述井筒内气相密度沿井筒深度分布的结果由下式表示：

其中，ρ_g为气相密度，p为气液固三相混合物的压力，T为井筒温度，Z为气相偏差系数，R为通用气体常数，M为甲烷摩尔质量。

根据本发明的一个实施例，在所述连续的子井段内根据所述压力模型和温度分布模型耦合迭代计算直至井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果包括：

耦合计算步骤，依据当前子井段上出口位置的压力和温度，液相体积流量参数和气相体积流量参数，根据所述压力模型和温度分布模型耦合计算，获得当前子井段下入口位置的气液固三相流互相耦合的压力和温度，以及当前子井段下入口位置的液相体积流量参数和气相体积流量参数；

迭代计算步骤，以当前子井段下入口位置的压力和温度做为下一子井段上出口位置的压力和温度，以当前子井段下入口位置的液相体积流量参数和气相体积流量参数做为下一子井段上出口位置的液相体积流量参数和气相体积流量参数；

重复执行耦合计算步骤和迭代计算步骤，直到井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

根据本发明的一个实施例，所述所述耦合计算步骤包括，

设定当前子井段下入口的预设温度，依据当前子井段上出口位置的液相体积流量参数和气相体积流量参数，根据所述气液固三相流压力模型计算当前子井段下入口位置的压力；由所述当前子井段下入口位置的压力和所述温度分布模型得到当前井段下入口的计算温度；

当所述预设温度和计算温度之间的差值小于预设的温度差阈值时，所述压力和预设温度做为当前井段下入口位置互相耦合的压力和温度，并计算当前子井段下入口位置的液相体积流量参数和气相体积流量参数；当所述预设温度和计算温度之间的差值大于预设的温度差阈值时，重新设定当前子井段下入口的预设温度。

本发明带来了以下有益效果：可以根据井口产液量和产气量，以及套压等物性参数预测任意井型的井筒内任意位置的流态，流速，压力，温度分布等物性参数；可以进行生产参数敏感性分析，依据不同的井底压力预测单井产量变化，然后依据配产需求进行参数调整。

本发明的其它特征和优点将在随后的说明书中阐述，并且部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

图1是煤层气井井筒流动状态示意图；

图2是根据本发明的一个实施例的煤层气井井筒流动动态预测方法流程图；

图3是根据本发明的一个实施例的耦合迭代计算的步骤流程图；

图4是根据本发明的另一个实施例的耦合迭代计算的步骤流程图；

图5是根据本发明实施例的方法获得的井筒压力沿井筒轴向的分布曲线；

图6是在工程现场得到的油套环空压力分布曲线；

图7是在工程现场得到的油套环空温度分布曲线。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

首先对本发明的应用环境进行说明。图1所示为油管产水套管产气生产方式下煤层气井井筒流动动状态示意图。煤层气储层初始情况下含水率通常较高，同时由于大规模压裂改造以及孔隙壁面摩擦导致在排采过程中通常伴随有煤粉颗粒析出，因此煤层气井筒油管101内为固液两相流，油套环空102内动液面以上为纯气段，动液面至泵吸入口103附近可以近似认为是静液柱段，泵吸入口103以下的套管104内为气液固三相流，其中固相煤粉颗粒为分散相。

由于井筒内流体的物性参数通常是耦合在一起的，即压力、温度数值会相互影响，须将压力和温度参数进行耦合求解。本发明的实施例提供一种综合压降计算模型来模拟井筒内气液固三相流的流动形态，即针对不同流动型态的特点，分别采用与之相对应的压降计算模型，从而克服以往采用单一模型计算压差造成的计算精度不高的弊端。本发明实施例提供的综合压降计算模型优点更接近于实际，提高计算精度。

以下结合具体实施例对本发明提供的方法和装置进行详细说明。

实施例一

本发明公开了一种煤层气井井筒流动动态预测方法，图2所示为流程图。以下结合图2对本实施例的方法作详细说明。

在步骤S201中，对井筒油管中的液体采样以获取井筒中的真实固体含量H_s；测量井口的液相体积流量和气相体积流量，以及井口温度和压力。

具体的，在油管产水套管产气排采方式下分别从油管取水样，固相为煤粉，煤粉浓度从水样中测出。若取出的水样中煤粉颗粒小，绝大部分能够悬浮于水中且沉淀极少时，可使用浊度计测量悬浮煤粉颗粒浓度。若煤粉颗粒较大造成沉淀时，可通过蒸馏法获取一定体积水样内干煤粉质量，之后换算成煤粉浓度，从而得到真实固体含量H_s。

可由设置在生产现场井口的流量计测量得到井口的液相体积流量和气相体积流量，由设置在油管口和油套环空口的压力计和温度计测量井口温度和压力。

在步骤S202中，确定井筒中气液固三相流的流动形态，并根据井筒中的真实固体含量H_s确定气液固三相流的建模参数。

其中，所述建模参数包括井筒内的真实液体含量H_l(θ)，真实气体含量H_g(θ)，沿程阻力系数λ。

表1

由于本发明的实施例引入真实固体含量H_s，固相密度ρ_s来表示固相煤粉的参数，需要对现有气液两相流分析方法中的液体含量H_l'(θ)进行修正；同时，气液固三相流沿程阻力系数λ也需在考虑煤粉颗粒的情况下进行修正。

具体的，本步骤提出了针对气液固三相流的流态判别准则，假设固相颗粒在井筒任意截面上的平均浓度保持不变，即真实固体含量H_s在井筒任意深度处保持不变，且固相和液相在垂直方向上无相对运动。应用贝格斯-布里尔(Beggs-Brill)方法，采用表1中的界限来划分各种流态，而各参数的计算方法如下所述。

其中，弗洛德准数：

$N_{\mathrm{Fr}}=\frac{v_m^2}{\mathrm{gD}}---\left(1\right)$

(1)式中g为重力加速度，D为管道直径，v_m为混合物平均流速。

无滑脱持液率：

$E_L=\frac{Q_l}{Q_l+Q_g}---\left(2\right)$

式中，Q_l为井筒中的液相体积流量；

Q_g为井筒中的气相体积流量。

L₁，L₂，L₃和L₄为四个流型区的分隔线，分区线的方程为：

$\begin{array}{c}{L}_1=316E_L^{0.302}\\ L_2=0.0009252{E_L}^{-2.4684}\\ L_3=0.10E_L^{-6.733}\\ L_4=0.5E_L^{-6.733}\end{array}$

根据Q_l、Q_g便可计算混合物平均流速和气相表观流速：

混合物平均流速： $v_m=\frac{Q_l+Q_g}{0.25\mathrm{\π}{D}^2}---\left(3\right)$

气相表观流速： $v_{\mathrm{sg}}=\frac{Q_g}{0.25\mathrm{\π}{D}^2}---\left(4\right)$

上式中，D为管道直径。

在用Beggs-Brill方法进行计算倾斜管流时，首先按水平管计算,然后进行倾斜角校正。

H'_l(θ)＝H_l(0)ψ(5)

(5)式中，H'_l(θ)为倾角为θ的气液两相流动的液体含量；Hl(0)为同样流动参数下,水平流动时的液体含量；ψ为倾斜校正系数。

$H_l\left(0\right)=\frac{a{E}_L^b}{N_{\mathrm{Fr}}^C}---\left(6\right)$

式中，a、b、c为取决于流型的常数(见表2)

表2

利用表2和上式计算出的H_l(0)必须满足H_l(0)≥E_L，否则,取H_l(0)＝E_L。因为E_L是无滑脱时的持液率，而H_l(0)为存在滑脱时的持液率，因此，H_l(0)的最小值是E_L。

实验结果表明，倾斜校正系数ψ不仅与倾斜角θ有关，而且与无滑脱持液率E_L、弗洛德数N_Fr及液相速度数N_vl有关。

根据实验结果回归的倾斜校正系数ψ的相关式如下:

$\mathrm{\ψ}=1+C[\sin \left(1.8\mathrm{\θ}\right)-\frac{1}{3}{\sin }^3\left(1.8\mathrm{\θ}\right)]---\left(7\right)$

对于垂直管：

ψ＝1+0.3C

系数C与无滑脱持液率E_L、弗洛德数N_Fr和液相速度数N_vl有关。

$N_{\mathrm{vl}}=V_{\mathrm{sl}}{\left(\frac{{\mathrm{\ρ}}_l}{\mathrm{g\σ}}\right)}^{1/4}---\left(8\right)$

式中，为液相表观流速；σ为液体表面张力；g为重力加速度。

C＝(1-E_L)ln[d(E_L)^e(N_vl)^f(N_Fr)^g1](9)

(9)式中的系数d、e、f和g₁由表3根据流型来确定。

表3

由于流体中还有煤粉固体颗粒存在，因此需修正气体和液体的真实含量，修正后的液体真实含量和气体真实含量分别为：

H_l(θ)＝(1-H_s)H'_l(θ)(10)

H_g(θ)＝(1-H_s)[1-H'_l(θ)](11)

以下步骤确定气液固三相流相流沿程阻力系数λ：

沿程阻力系数

λ＝λ′·e^s(12)式中

λ′为无滑脱的沿程阻力系数，无因次；s为指数。

上式中的λ′可由下式计算

${\mathrm{\λ}}^{\′}{\left[21g\left(\frac{R_e^{\′}}{4.52231g{R}_e^{\′}-3.8215}\right)\right]}^{-2}---\left(13\right)$

含有固体小颗粒的液相粘度为

μ_ls＝μ_l(1+2.5H_s)(14)

结合上式可得

$R_e^{\′}=\frac{\mathrm{Dv}[{\mathrm{\ρ}}_l{H}_l+{\mathrm{\ρ}}_g{H}_g+{\mathrm{\ρ}}_s{H}_s]}{{\mathrm{\μ}}_{\mathrm{ls}}(H_l+H_s)+{\mathrm{\μ}}_g{H}_g}---\left(15\right)$

式中

R_e′为无滑脱的雷诺数；μ_l、μ_g分别为液相、气相的粘度，单位：Pa·s。而

$s=\frac{\ln Y}{-0.0523+3.182\ln Y-0.8725{\left(\ln Y\right)}^2+0.01853{\left(\ln Y\right)}^4}---\left(16\right)$

其中

$Y=\frac{E_L}{{\left[H_l\left(\mathrm{\θ}\right)\right]}^2}---\left(17\right)$

需要指出，当1<Y<1.2时，依据现有理论应使用下式求s

s＝ln(2.2Y-1.2)(18)

在步骤S203中，依据所述真实固体含量H_s以及所述建模参数，结合动液面位置和井身结构数据，基于贝格斯-布里尔方法建立井筒内气液固三相流压力模型。

所述气液固三相流压力模型由下式表示：

$\frac{\mathrm{dp}}{\mathrm{dz}}-\frac{[{\mathrm{\&rho;}}_l{H}_l\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_g{H}_g\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_s{H}_s]g\sin \mathrm{\&theta;}+\frac{\mathrm{\&lambda;}{\mathrm{Gv}}_m}{2\mathrm{DA}}}{1-\frac{[{\mathrm{\&rho;}}_l{H}_l\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_g{H}_g\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_s{H}_s]v_m{v}_{\mathrm{sg}}}{p}}---\left(19\right)$

上式中：

p为混合物的绝对压力，单位：Pa；z为沿井筒轴向流动的距离，单位：m；ρ_l为液相密度，单位：kg/m³；ρ_g为气相密度，单位：kg/m³；ρ_s为固相密度，即煤粉的密度，单位：kg/m³；H_l(θ)为真实液体含量，单位：m³/m³；H_g(θ)为真实气体含量，单位：m³/m³；H_s为真实固体含量，单位：m³/m³；g为重力加速度，单位：m/s²；θ为井筒管道与水平方向的夹角，单位：°；λ为气液固三相流沿程阻力系数，无因次；G为混合物的质量流量，单位：kg/s；v_m为子井段上出口位置的混合物的平均流速，单位：m/s；v_sg为子井段上出口位置的气相表观流速，单位：m/s；D为管道直径，单位：m；A为管道截面积，单位：m²。为简洁起见，本实施例中的混合物均表示井筒内气液固三相混合物。

以下给出混合物的质量流量G的计算过程。

现场生产过程中井口油管产出液相和煤粉，油套环空产出煤层气体。现场通过流计量测量出液体和气体的井口体积流量Q_l0、Q_g0，以及井口压力和温度p0、T0。已知液相密度ρ_l(其中ρ^l为常数，将在下文中解释),同时依据公式(20)计算出井口气相密度ρ_g0后，便可计算出液相和气相的质量流量：

G_l0＝ρ_lQ_l0，G_g0＝ρ_g0Q_g0

根据井口油管口测得煤粉含量H_s可计算得煤粉质量流量：

G_s0＝ρ_sQ_l0H_s

由于生产过程中井筒内流动可在瞬时达到稳态过程，压力、温度和质量流量参数与时间无关，即在任一位置，气液固三相的质量流量均分别为G_g0、G_l0和G_s0。因此三相混合物密度为G＝G_l0+G_g0+G_s0。

进一步，煤层气井中液相为水，水的密度受温度、压力影响不大，因此本发明中ρ_l取水在标准状态下的密度1000kg/m³。

在油管产水套管产气排采方式下分别从套管取煤层气样。使用气相色谱仪测量煤层气体组分和浓度，由于煤层气中甲烷含量最高，通常在90％以上，因此这里近似认为气相为甲烷。依据气体状态方程，气相密度为：

${\mathrm{\&rho;}}_g=\frac{\mathrm{pM}}{\mathrm{ZRT}}---\left(20\right)$

式中，T为混合物温度，Z为气相偏差系数，R为通用气体常数8.314Pa.m³/(mol.K)，M为甲烷摩尔质量16g/mol。本发明中使用Hall-Yarborough方法求解Z：

$Z=\left[\frac{0.06125p_r{T}_r^{-1}}{y}\right]\exp [-1.2{(1-T_r^{-1})}^2]---\left(21\right)$

y为下列方程(4)的解

$\begin{array}{c}-0.06125p_r{T}_r^{-1}\exp [-1.2{(1-T_r^{-1})}^2]+\frac{y+y^2+y^3-y^4}{{(1-y)}^3}-({14.76T}_r^{-1}-{9.76T}_r^{-2}+4.58T_r^{-3})y^2\\ +({90.7T}_r^{-1}-{242.2T}_r^{-2}+42.4T_r^{-3})y^{(2.18+{2.82T}_r^{-1})}=0\end{array}---\left(22\right)$

公式(3)、(4)、(21)、(22)中，p_r＝p/p_c为对比压力，T_r＝T/T_c为对比温度，p为混合物压力，T为混合物温度。对于甲烷，p_c＝4.6408MPa，T_c＝190.67K。

根据井型不同，煤层气井各段θ值不同。对于直井和斜井，各段θ值均可依据钻井数据直接获取；对于水平井，可依据井眼轨迹数据，使用自然参数法确定各段θ值。

上式(19)给出气液固三相流动压力模型。结合图1，对于油套环空102内动液面以上的纯气段以及油管101内的液固两相流均可以看做三相流的特殊情况。例如：纯气段计算可以通过假设液相和固相含量为零并继续使用上述模型进行。即：对于油套环空内动液面以上的纯气段，H_l(θ)＝0,H_s＝0；对于油管内的固液两相流，H_g(θ)＝0。该处理方法在上述步骤S202中流型判别等模型中同样适用。

至此为止，建立了煤层气井油管产水套管产气时可压缩多组分流体井筒气液固三相流动压降计算模型，给出了流态判别准则及相应阻力系数的计算方法，解决了任意井型煤层气井井筒压力计算问题。

在步骤S204中，依据所述气液固三相流压力模型，所述真实固体含量H_s以及所述建模参数，结合井身传热参数，根据能量守恒性质建立子井段内温度分布模型。

具体的，煤层气井井筒内能量守恒方程：

$\frac{\mathrm{dh}}{\mathrm{dz}}=\frac{\mathrm{dq}}{\mathrm{dz}}-\frac{v_m{\mathrm{dv}}_m}{\mathrm{dz}}-g\sin \mathrm{\&theta;}-\frac{{{\mathrm{\&lambda;v}}_m}^2}{2D}---\left(23\right)$

其中，h为比焓，q为径向热流量，D为管道直径，λ为步骤S202中确定的气液固三相流相流沿程阻力系数，v_m为步骤S202中得到的混合物平均流速。

假设煤层气的流动过程为等焓过程，得到：

dh＝C_pmdT-C_pmC_Jmdp(24)

其中，C_pm为混合物的平均定压比热容，C_Jm为混合物的焦耳—汤姆逊数。

代入能量方程得到井筒内温度分布模型：

$C_{\mathrm{pm}}\frac{\mathrm{dT}}{\mathrm{dz}}-C_{\mathrm{pm}}{C}_{\mathrm{Jm}}\frac{\mathrm{dp}}{\mathrm{dz}}+\frac{v_{m}d{v}_m}{\mathrm{dz}}+g\sin \mathrm{\&theta;}+\frac{\mathrm{\&lambda;}{v}_m^2}{2D}=-\frac{\mathrm{dq}}{\mathrm{dz}}---\left(25\right)$

由于生产过程中井筒内流动和传热可在瞬时达到稳态过程，流动参数与传热参数与时间无关，因此井筒传热是稳定的,于是煤层气井井筒内热传导方程：

$\frac{\mathrm{dq}}{\mathrm{dz}}=-\frac{2{\mathrm{\&pi;r}}_{\mathrm{to}}{U}_{\mathrm{to}}}{G}(T-T_h)---\left(26\right)$

式中：T_h为井壁温度；G为气液固三相混合物的质量流量；U_to为井筒总传热系数。

应用Ramey推荐的无因次时间函数f(t_D)，上式可表示为

$\frac{\mathrm{dq}}{\mathrm{dz}}=-\frac{2\mathrm{\&pi;}{K}_e}{\mathrm{Gf}\left(t_D\right)}(T_h-T_e)---\left(27\right)$

式中K_e、T_e分别为地层传热系数、地层初始温度；

用Hasan-Kabir公式计算f(t_D)，

$f\left(t_D\right)=\left\{\begin{array}{cc}(0.5\ln t_D+0.4063)(1+0.6/t_D)& (t_D>1.5)\\ 1.1281\sqrt{t_D}(1-0.3\sqrt{t_D})& ({10}^{-10}\&le;t_D\&le;1.5)\end{array}\right.---\left(28\right)$

其中，α为地层热扩散系数；t为油井生产时间；r_wb为井眼半径。

联立上式消除井壁温度T_h，可得热流梯度方程，

$\frac{\mathrm{dq}}{\mathrm{dz}}=-\frac{2\mathrm{\&pi;}{r}_{\mathrm{to}}{U}_{\mathrm{to}}{K}_e}{W_m[r_{\mathrm{to}}{U}_{\mathrm{to}}f\left(t_D\right)+K_e]}(T_f-T_e)---\left(29\right)$

井筒总传热系数的计算公式为

$U_{\mathrm{to}}={[\frac{r_{\mathrm{to}}}{r_{\mathrm{ti}}{h}_f}+\frac{r_{\mathrm{to}}\ln (r_{\mathrm{to}}/r_{\mathrm{co}})}{K_t}+\frac{1}{h_r+h_c}+\frac{r_{\mathrm{to}}\ln (r_{\mathrm{co}}/r_{\mathrm{ci}})}{K_c}+\frac{r_{\mathrm{to}}\ln (r_{\mathrm{wb}}/r_{\mathrm{co}})}{K_{\mathrm{cem}}}]}^{-1}---\left(30\right)$

式中：

r_ci、r_co为套管内、外半径，单位：m；

r_ti、r_to为油管内、外半径，单位：m；

h_f为管柱内流体与油管之间的传热系数，单位：W/(m².K)；

h_c、h_r为环空气体对流和辐射传热系数，单位：W/(m².K)；

K_cem为水泥环的导热系数，单位：W/(m.K)；

K_t、K_c为油管、套管水泥环的导热系数,单位：W/(m.K)。

至此为止，已建立井筒内温度分布模型。

在步骤S205中，从井口开始对井筒依次划分为若干连续的子井段，将测量得到的井口的液相体积流量和气相体积流量，以及井口温度和压力做为初始值，在所述连续的子井段内根据所述压力模型和温度分布模型耦合迭代计算直至井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

优选的，还包括在步骤S206(图中未示出)中：根据所述压力和温度沿井筒深度的分布结果以及气体状态方程获得井筒内气相密度沿井筒深度的分布结果。参考公式(20)，所述井筒内气相密度沿井筒深度分布的结果由下式表示：

${\mathrm{\&rho;}}_g=\frac{\mathrm{pM}}{\mathrm{ZRT}},$

其中，ρ_g为气相密度，p为气液固三相混合物的压力，T为井筒温度，Z为气相偏差系数，R为通用气体常数，M为甲烷摩尔质量。

以下结合图3对在连续的子井段内根据所述压力模型和温度分布模型迭代计算直至井底的步骤流程作详细说明。

在耦合计算步骤S301中，由当前子井段上出口位置的压力和温度，当前子井段上出口位置的液相体积流量参数和气相体积流量参数，根据所述压力模型和温度分布模型耦合计算，获得当前子井段下入口位置的气液固三相流互相耦合的压力和温度，以及当前子井段下入口位置的液相体积流量参数和气相体积流量参数。

具体的，步骤S301包括以下子步骤：

子步骤S3011，设定当前子井段下入口的预设温度，依据当前子井段上出口位置的液相体积流量参数和气相体积流量参数，根据所述压力模型计算当前子井段下入口位置的压力；由所述当前子井段下入口位置的压力和所述温度分布模型得到当前井段下入口的计算温度；

子步骤S3012，比较所述预设温度和计算温度，判断所述预设温度和计算温度之间的差值是否小于预设的温度差阈值，若是，执行子步骤S3013；若否，执行子步骤S3011，重新设定当前子井段下入口的预设温度；

子步骤S3013，将所述压力和预设温度做为当前井段下入口位置互相耦合的压力和温度，并计算当前子井段下入口位置的液相体积流量参数和气相体积流量参数。这里依据质量守恒，即下入口位置与上出口位置三相流的质量流量不变，在此基础上依据下入口处的压力、温度计算结果可得出气液的密度，质量流量除以密度得到体积流量。具体的，下入口处液相密度为ρ_l，质量流量为G_lin＝G_l0，气相密度为ρ_gin，质量流量为G_gin＝G_g0，则下入口处液相和气相体积流量分别为：

$Q_{\mathrm{lin}}=\frac{G_{\mathrm{lin}}}{{\mathrm{\&rho;}}_l},Q_{\mathrm{gin}}=\frac{G_{\mathrm{gin}}}{{\mathrm{\&rho;}}_{\mathrm{gin}}}.$

在迭代计算步骤S302中，以当前子井段下入口位置的压力和温度做为下一子井段上出口位置的压力和温度，以当前子井段下入口位置的液相体积流量参数和气相体积流量参数做为下一子井段上出口位置的液相体积流量参数和气相体积流量参数。

反复执行耦合计算步骤S301和迭代计算步骤S302，直到井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

优选的，步骤S205可以由图4所示的一系列子步骤完成。包括：

子步骤S2051，开始步骤，设置温度差阈值ε；从井口开始对井筒依次划分为N个连续的子井段，从井口开始依次记为第1子井段，第2子井段，…，第N子井段；第1子井段上出口位置的液相体积流量Q_l和气相体积流量Q_g由设置在井口的流量计测量得到；第1子井段上出口位置的温度和压力由设置在油管口和油套环空口的压力计和温度计测量；

子步骤S2052，设定第i子井段下入口的预设温度其中i为大于等于1的整数，i初始值为1；获取第i子井段上出口位置的液相体积流量Q_l和气相体积流量Q_g；

子步骤S2053，根据所述压力模型获得第i子井段下入口位置的压力根据和所述温度分布模型得到第i子井段下入口位置的计算温度

子步骤S2054，判断预设温度与计算温度之间的温度差是否小于温度差阈值，即判断是否若是，则执行子步骤S2055；若否，则执行子步骤S2052，重新设定预设温度

子步骤S2055，将压力和预设温度分别做为第i子井段下入口位置的互相耦合的压力P_in和温度T_in，即计算第i子井段下入口位置的液相体积流量Q_lin和气相体积流量Q_gin；

已知下入口位置压力P_in和温度T_in，依据方程(20)可算出下入口位置气相密度ρ_gin。下入口处液相密度为ρ_l，质量流量为G_lin＝G_l0，气相密度为ρ_gin，质量流量为G_gin＝G_g0，则下入口处液相和气相体积流量分别为：

$Q_{\mathrm{lin}}=\frac{G_{\mathrm{lin}}}{{\mathrm{\&rho;}}_l},Q_{\mathrm{gin}}=\frac{G_{\mathrm{gin}}}{{\mathrm{\&rho;}}_{\mathrm{gin}}};$

将第i子井段下入口位置的液相体积流量Q_lin和气相体积流量Q_gin做为第i+1子井段上出口位置的液相体积流量Q_l和气相体积流量Q_g；

子步骤S2056，i数值加1；

子步骤S2057，判断是否到井底，若i≤N，表示未到井底，执行子步骤子步骤S2052；若i>N表示已到井底，结束，获得井筒内气液固三相流的压力和温度沿井筒深度的分布。

以下对井口温度和井口压力作为初始值的情况作详细说明：

参照图1，对于油管内的固液两相流，在油管101的油管口测量得到的油管口的温度和压力作为初始值，用于计算油管101内从泵吸入口103到油管口的固液两相流的压力和温度分布；

对于油套环空102内动液面以上的纯气段，在油套环空102的出口位置测量得到的温度和套压作为初始值，用于计算油套环空102内从动液面105到油套环空口的纯气段的压力和温度分布；

对于套管104内泵吸入口103位置的气液固三相流，其压力为动液面105至泵吸入口103附近的静液柱的压力与油套环空102内纯气段在动液面105位置的压力之和。

本发明的实施例可提供一种结合综合压降计算模型和流体热力性质的气液固三相流动态预测方法，从而确定压力、温度耦合的情况下煤层气井井筒物性参数，定量计算煤层气井井筒内压力、温度、密度等参数沿井筒轴向的分布值。

根据本实施例提供的方法可以根据井口产液量和产气量，以及套压等物性参数预测任意井型的井筒内任意位置的流态，流速，压力，温度分布等物性参数。图5所示为在不同套压的情况下井筒压力沿井筒轴向的分布，具体而言，起点是油套环空井口，一直沿油套环空向下到达动液面直到井底。

因此，依据本实施例提供的方法可以进行生产参数敏感性分析，如在不同管柱尺寸、套压等情况下的压力沿井筒分布，并依据不同的井底压力预测单井产量变化，然后依据配产需求进行参数调整。

本发明的实施例提供的方法在生产中可取得显著效果，主要体现在：

(1)生产现场通过回声仪测量动液面高度，根据现场产液量和产气量的变化，精确预测每口井井底压力的变化情况；根据某一区块多口井的井底压力变化情况，结合油藏软件进行历史拟合，最终确定该区块的实际可采储量，即整个区块的最终产气量；

(2)煤储层应力敏感性较强，排采过程中应尽量维持井底压力缓慢下降，避免压力在短时间内下降过快导致储层伤害，影响累计产量；根据本发明实施例提供的方法，可实时掌握井底压力值，并采取措施保持井底压力稳定降低；例如，若发现压力下降过快，可通过使用更小的油嘴降低产液量和产气量的方法增大井口压力，降低井底压力下降速度。

实施例二

本实施例为在工程现场应用本发明的方法的一个实例。该井为一口直井，井口测量煤粉平均体积浓度为5％，即真实固体含量H_s＝5％；井深为1399.49m；管道直径D＝139.7-7.72＝131.98mm；套压为1MPa；动液面高度806.32m。

工作现场井身数据表如表4所示。

表4

生产数据如表5所示。

表5

井身传热参数如表6所示。

表6

项目	值	单位	项目	值	单位
						环空气体传热系数	0.3	W/m2.K	管壁传热系数	1.73	W/m2.K
地层导热系数	1.717		环空气体辐射系数	0.22	W/m2.K
						气相导热系数	0.3		水泥环导热系数	0.57
液相定压比热	4.2		地层扩散系数	0.75	km2/s

根据实施例一提供的方法，可得到如图6所示的井筒内油套环空压力分布，以及如图7所示的油套环空温度分布。使用本发明计算井底压力4.16MPa，与实测井底压力4.25MPa仅相差2.21％。可见本发明提供的方法可以对煤层气井井筒内气液固三相流的压力分布进行精确分析。

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种煤层气井井筒流动动态预测方法，其特征在于，包括：

对井筒油管中的液体采样以获取井筒中的真实固体含量H_s；测量井口的液相体积流量和气相体积流量，以及井口温度和井口压力；

判别井筒中气液固三相流的流动形态，并根据井筒中的真实固体含量H_s确定气液固三相流的建模参数；

依据所述真实固体含量H_s以及所述建模参数，结合动液面位置和井身结构数据，基于贝格斯-布里尔方法建立井筒内的气液固三相流压力模型；

依据所述气液固三相流压力模型，所述真实固体含量H_s以及所述建模参数，结合井身传热参数，根据能量守恒性质建立井筒内的温度分布模型；

从井口开始对井筒依次划分为若干连续的子井段，将测量得到的井口的液相体积流量和气相体积流量，以及井口温度和井口压力做为初始值，在所述连续的子井段内根据所述压力模型和温度分布模型耦合迭代计算直至井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

2.根据权利要求1所述的方法，其特征在于，还包括：

根据所述压力和温度沿井筒深度的分布结果以及气体状态方程获得井筒内气相密度沿井筒深度的分布结果。

3.根据权利要求2所述的方法，其特征在于，所述判别井筒中气液固三相流的流动形态包括：

由井筒中的液相体积流量参数和气相体积流量参数确定气液固混合物的平均流速 $v_m=\frac{Q_l+Q_g}{0.25\mathrm{\&pi;}{D}^2}$ 和无滑脱持液率 $E_L=\frac{Q_l}{Q_l+Q_g};$

根据弗洛德准数和无滑脱持液率应用贝格斯-布里尔方法划分井筒中气液固三相流的流动形态；

其中，Q_l为井筒中的液相体积流量参数，Q_g为井筒中的气相体积流量参数，g为重力加速度，D为管道直径，v_m为气液固三相混合物的平均流速。

4.根据权利要求3所述的方法，其特征在于，所述根据真实固体含量H_s确定气液固三相流的建模参数包括：

根据真实固体含量H_s和井筒中气液固三相流的流动形态的流型参数计算真实液体含量H_l(θ)＝(1-H_s)H'_l(θ)，真实气体含量H_g(θ)＝(1-H_s)[1-H'_l(θ)]，

其中，H'_l(θ)＝H_l(0)ψ，H'_l(θ)为倾角为θ的气液两相流动的液体含量，H_l(0)为同样流型参数下水平流的液体含量，ψ为倾斜校正系数，θ为井筒管道与水平方向的夹角。

5.根据权利要求4所述的方法，其特征在于，所述根据真实固体含量H_s确定气液固三相流的建模参数还包括：

根据真实固体含量H_s，真实液体含量H_l(θ)和无滑脱持液率E_L计算气液固三相流沿程阻力系数λ＝λ′·e^s；其中，

无滑脱的沿程阻力系数 ${\mathrm{\&lambda;}}^{\&prime;}={\left[21g\left(\frac{R_e^{\&prime;}}{4.52231g{R}_e^{\&prime;}-3.8215}\right)\right]}^{-2};$

R_e′为无滑脱的雷诺数，由含有固相颗粒的液相粘度μ_ls＝μ_l(1+2.5H_s)确定，μ_l为液相粘度，H_s为真实固体含量；

指数 $s=\frac{\ln Y}{-0.0523+3.182\ln Y-0.8725{\left(\ln Y\right)}^2+0.01853{\left(\ln Y\right)}^4},$

其中，E_L为无滑脱持液率，H_l(θ)为真实液体含量，θ为井筒管道与水平方向的夹角。

6.根据权利要求5所述的方法，其特征在于，

所述气液固三相流压力模型由下式表示：

$\frac{\mathrm{dp}}{\mathrm{dz}}-\frac{[{\mathrm{\&rho;}}_l{H}_l\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_g{H}_g\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_s{H}_s]g\sin \mathrm{\&theta;}+\frac{\mathrm{\&lambda;}{\mathrm{Gv}}_m}{2\mathrm{DA}}}{1-\frac{[{\mathrm{\&rho;}}_l{H}_l\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_g{H}_g\left(\mathrm{\&theta;}\right)+{\mathrm{\&rho;}}_s{H}_s]v_m{v}_{\mathrm{sg}}}{p}},$

其中，ρ_l为液相密度，ρ_g为气相密度，ρ_s为固相密度，p为气液固三相混合物的压力，z为沿井筒轴向流动的距离，g为重力加速度，G为气液固三相混合物的质量流量，A为管道横截面积，D为管道直径，v_m为气液固三相混合物的平均流速，v_sg为气相表观流速；H_s为真实固体含量，H_l(θ)为真实液体含量，H_g(θ)为真实气体含量，θ为井筒管道与水平方向的夹角；λ为沿程阻力系数。

7.根据权利要求6所述的方法，其特征在于，

所述温度分布模型由下式表示：

$C_{\mathrm{pm}}\frac{\mathrm{dT}}{\mathrm{dz}}-C_{\mathrm{pm}}{C}_{\mathrm{Jm}}\frac{\mathrm{dp}}{\mathrm{dz}}+\frac{v_{m}d{v}_m}{\mathrm{dz}}+g\sin \mathrm{\&theta;}+\frac{\mathrm{\&lambda;}{v}_m^2}{2D}=-\frac{\mathrm{dq}}{\mathrm{dz}}$

其中，由所述气液固三相流压力模型确定，T为井筒温度，C_pm为气液固三相混合物的平均定压比热容，C_Jm为气液固三相混合物的焦耳—汤姆逊数，q为径向热流量，p为气液固三相混合物的压力，v_m为气液固三相混合物的平均流速，λ为气液固三相流沿程阻力系数，θ为井筒管道与水平方向的夹角，D为管道直径，g为重力加速度，z为沿井筒轴向流动的距离。

8.根据权利要求7所述的方法，其特征在于，所述井筒内气相密度沿井筒深度分布的结果由下式表示：

其中，ρ_g为气相密度，p为气液固三相混合物的压力，T为井筒温度，Z为气相偏差系数，R为通用气体常数，M为甲烷摩尔质量。

9.根据权利要求8所述的方法，其特征在于，在所述连续的子井段内根据所述压力模型和温度分布模型耦合迭代计算直至井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果包括：

耦合计算步骤，依据当前子井段上出口位置的压力和温度，液相体积流量参数和气相体积流量参数，根据所述压力模型和温度分布模型耦合计算，获得当前子井段下入口位置的气液固三相流互相耦合的压力和温度，以及当前子井段下入口位置的液相体积流量参数和气相体积流量参数；

迭代计算步骤，以当前子井段下入口位置的压力和温度做为下一子井段上出口位置的压力和温度，以当前子井段下入口位置的液相体积流量参数和气相体积流量参数做为下一子井段上出口位置的液相体积流量参数和气相体积流量参数；

重复执行耦合计算步骤和迭代计算步骤，直到井底，获得井筒内气液固三相流的压力和温度沿井筒深度的分布结果。

10.根据权利要求9所述的方法，其特征在于，所述耦合计算步骤包括，

设定当前子井段下入口的预设温度，依据当前子井段上出口位置的液相体积流量参数和气相体积流量参数，根据所述气液固三相流压力模型计算当前子井段下入口位置的压力；由所述当前子井段下入口位置的压力和所述温度分布模型得到当前井段下入口的计算温度；

当所述预设温度和计算温度之间的差值小于预设的温度差阈值时，所述压力和预设温度做为当前井段下入口位置互相耦合的压力和温度，并计算当前子井段下入口位置的液相体积流量参数和气相体积流量参数；当所述预设温度和计算温度之间的差值大于预设的温度差阈值时，重新设定当前子井段下入口的预设温度。