产气井温度压力分布的数值模拟方法

摘要

本发明涉及油气藏开发工程管理技术领域，其公开了一种产气井温度压力分布的数值模拟方法，对产气井的温度压力参数分布情况进行准确模拟，从而为生产提供指导依据。该方法包括以下步骤：A、建立温度压力耦合微分方程组模型；B、估计微分方程组相关参数；C、使用差分法进行算法设计；D、基于C#开发平台构建模拟软件；E、在模拟软件中输入初始数据，获得温度压力分布模拟结果。本发明能够对产气井的温度和压力分布进行精确预测，极大促进油气开采设备设计水平的提高，从而有利于油气藏开发。

说明书

技术领域

本发明涉及油气藏开发工程管理技术领域，具体为产气井温度压力分布的数值模拟方法， 涉及产气井的温度压力耦合机理分析，数学模型建立，数值模拟方法设计等。

背景技术

生产气井通常伴随着与高温、高流速、高压力，这是由于摩擦、井筒形变、热传导等原 因。随着生产环境的演化，包括深水和高温度压力条件，在石油工程设施的设计、水合物预 防优化、生产井的动态分析，有必要对温度压力进行准的预测。

当液体从一个区域最初产生时，其温度在井底可能认为一样的。这对气体是不成立的。 如果Joule-Thompson效应被恰当考虑，入口气体温度可用地层温度进行估计。因此，井筒底 部的温度，能够被可靠地估计。然而，随着流体的上升，它的温度显著高于周围地层的温度， 这是因为地层温度随着深度减少而下降。

当地层与流体有温度差，将发生热传导现象。在任何深度，地层温度不仅与辐射距离有 关，也与生产时间有关。当达到稳定流时，湍流保证了在一定的深度，流体的温度是常数。 于是，流体中的热量损失随着时间而减少，并取决于井筒红热液与周围地层的不同的热阻。

完整的系统由流体，含有低压空气的环形间隙、套管、座封和地层，如图1所示：

导管内径是r_ti，外径是r_t0，套管的内径是r_ci，外径是r_c0，热量通过环圈中的空气进行传导。 辐射与对流也发生。当管体被加热，辐射能的传导速率取决于管体的温度。导管与套管之间 的辐射能的传导，取决于界面发出和吸收热量的特征。很多情况下，导管与孔之间的缝隙是 被封堵的。由于水泥的传导性可能低于周边的地层，计算由井筒地步开始，逐步向上。将气 井的进气口作为坐标轴的交点，井筒向下为坐标轴的正方向。图2显示了井筒微元的压减分 析，P是流体压力，v是气体速度，l是深度，dv是dl上的速度增量，dp是dl上的压力增 量，θ是井筒的倾斜角。

发明内容

本发明所要解决的技术问题是：提出一种产气井温度压力分布的数值模拟方法，对产气 井的温度压力参数分布情况进行准确模拟，从而为生产提供指导依据。

本发明解决上述技术问题所采用的方案是：产气井温度压力分布的数值模拟方法，包括 以下步骤：

A、建立温度压力耦合微分方程组模型；

B、估计微分方程组相关参数；

C、使用差分法进行算法设计；

D、基于C#开发平台构建模拟软件；

E、在模拟软件中输入初始数据，获得温度压力分布模拟结果。

进一步的，步骤A中，所述建立温度压力耦合微分方程组模型具体包括：

基于假设的模型构建条件分别获得能量守恒方程、质量守恒方程、动量守恒方程：

其中，能量守恒方程：

$q+\frac{d}{dl}(h+v^2/2-gl\sin \theta )=0---\left(1\right)$

质量守恒方程：

$\rho \frac{dv}{dl}+v\frac{d\rho }{dl}=0---\left(2\right)$

动量守恒方程：

$-\frac{1}{\rho }\frac{dP}{dl}+f\frac{v^2}{2r_{ti}}+gsin\theta =\frac{vdv}{dl}---\left(3\right)$

在微元dl上，从流体到座封-地层界面传导的辐射热为

$q=\frac{(-2\pi )r_{ti}{U}_{ti}}{w}(T-T_h)---\left(4\right)$

从座封-地层界面到周围地层传导的辐射热为：

$q=-\frac{2\pi ke}{wf\left(t_D\right)}(T_h-T_e)---\left(5\right)$

结合式(4)和式(5)得到流体与周围地层热传导的微分方程：

$q=\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\right](T-T_e)---\left(6\right)$

将(1)和(6)整合，得到如下常微分方程：

$\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}\left(t_D\right)+ke}\right](T-T_e)+\frac{d}{dl}(h+v^2/2-gl\sin \theta )=0---\left(7\right)$

当流体在井筒中流动的时候，由于管径变化很小，Joule–Thomson系数可以忽略不计，因 此dh＝C_pdT，对于式(7)整体的焓变化量为：

$\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\right](T-T_e)+C_p\frac{dT}{dl}+v\frac{dv}{dl}-g\sin \theta =0---\left(8\right)$

由式(1)、(2)、(3)、(8)可得天然气井中的流体压力、温度、密度和速率的微 分组方程：

$\{\begin{array}{c}\frac{d\rho }{dl}=\frac{\frac{-8.314z\rho }{C_{p}M}[\frac{\left(2\pi \right)r_{ti}{U}_{ti}ke}{w[r_{ti}{U}_{ti}f\left(t_D\right)+ke]}(T-T_e)-g\cos \theta ]-f\frac{v^2\rho }{2r_{ti}}-g\rho \cos \theta }{v^2-(\frac{8.314{\mathrm{zv}}^2}{C_{p}M}+\frac{8.314Tz}{M})}\\ \frac{dv}{dl}=-\frac{v}{\rho }\frac{d\rho }{dz}\\ \frac{dP}{dl}=\rho g\cos \theta +f\frac{v^2\rho }{2r_{ti}}+v^2\frac{d\rho }{dl}\\ \frac{dT}{dl}=\frac{[\frac{v^2}{\rho }\frac{d\rho }{dl}+g\cos \theta -\frac{\left(2\pi \right)r_{ti}{U}_{ti}ke}{w[r_{ti}{U}_{ti}f\left(t_D\right)+ke]}(T-Te)]}{C_p}\end{array}---\left(9\right)$

用y_i(i＝1,2,3,4)代替ρ，v，P，T，则方程组可简化为

$\frac{{\mathrm{dy}}_i}{dl}=F_i(y_1,y_2,y_3,y_4),(i=1,2,3,4)---\left(10\right).$

进一步的，所述假设的模型构建条件包括：

天然气的流动是单向稳定的；井筒的传热是稳态的；地层传热是不稳定的；导管和套管 是同心的。

进一步的，步骤B中，估计微分方程组相关参数具体包括：

B1.热传导系数估计：

$U_{ti}^{-1}=\frac{1}{hc+hr}+\frac{r_{ti}\ln \left(\frac{r_{cem}}{r_{co}}\right)}{k_{cem}}+\frac{r_{ti}\ln \left(\frac{r_{ci}}{r_{to}}\right)}{k_{ang}}---\left(11\right)$

B2.计算定压比热C_P：

C_p＝1243+3.14T+7.931×10^-4T²-6.881×10^-7T³(12)

B3.地层瞬态热传导方程：

$\left(\begin{array}{cc}f\left(t_D\right)=1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}& (t_D\le 1.5)\\ f\left(t_D\right)=(1+\frac{0.6}{t_D})[0.4063+0.5\ln \left(t_D\right)]& (t_D>1.5)\end{array}\right)---\left(13\right)$

$t_D=\frac{t\alpha }{r_{wb}^2}$

B4.计算摩擦系数f：

$\frac{1}{\sqrt{f}}=1.14-2\lg (\frac{0.00001524}{r_{ti}}+\frac{21.25}{{\mathrm{Re}}^{0.9}})---\left(14\right).$

进一步的，步骤C中，所述使用差分法进行算法设计包括对温度压力耦合微分方程组模 型的求解：

将产气井划分为n个区间，j(j＝1,2,…,n)为节点，n为井底，步长为h；

温度压力耦合微分方程组的边界条件为：

$P\left[n\right]=P_0,T\left[n\right]=T_0,\rho \left[n\right]=0.000001\times 3484.48{\gamma }_g\frac{P_0}{{\mathrm{zT}}_0},v\left[n\right]=\frac{101000\times 300000T_0}{293\times 86400P_{0}A}---\left(15\right)$

为了计算各节点的相关量，使用如下的迭代方程：

$y_i[j-1]=y_i\left[j\right]-\frac{h}{6}\left(a\right[i]+2b[i]+2c[i]+d[i\left]\right),(i=1,2,3,4;j=1,2,3,4...n)---\left(16\right)$

为使用以上方程，需要估计a[i],b[i],c[i]和d[i]:

$\left(\begin{array}{c}a\left[i\right]=F_i\left(y_1\right[n],y_2[n],y_3[n],y_4[n\left]\right)\\ b\left[i\right]=F_i\left(y_1\right[n]+\frac{h}{2}a[1],y_2[n]+\frac{h}{2}a[2],y_3[n]+\frac{h}{2}a[3],y_4[n]+\frac{h}{2}a[4\left]\right)\\ c\left[i\right]=F_i\left(y_1\right[n]+\frac{h}{2}b[1],y_2[n]+\frac{h}{2}b[2],y_3[n]+\frac{h}{2}b[3],y_4[n]+\frac{h}{2}b[4\left]\right)\\ d\left[i\right]=F_i\left(y_1\right[n]+hc[1],y_2[n]+hc[2],y_3[n]+hc[3],y_4[n]+hc[4\left]\right)\end{array}\right)---\left(17\right).$

进一步的，步骤D中，所述模拟软件包括用户界面模块、输出输出模块、算法模块和图 表模块；

所述用户界面模块用于提供气井的初始信息，这些信息可划分为静态数据、流体数据和 生产数据：

静态数据：描述静态下气井特征；

流体数据：描述流体的状态；

生产数据：生产阶段给定的参数；

所述输入输出模块用于从存储器里读入、保存数据；

所述算法模块用于根据设计的算法进行数学计算；

所述图表模块用于将模拟的数据结果进行图表化展示。

本发明的有益效果是：对产气井的温度和压力分布进行精确预测，极大促进油气开采设 备设计水平的提高，从而有利于油气藏开发。

附图说明

图1为井段结构示意图；

图2为导管内压减分析示意图；

图3为本发明产气井温度压力分布的数值模拟方法流程图；

图4为模拟软件结构示意图；

图5为程序流程示意图；

图6为模拟出来的压力分布曲线图；

图7为模拟出来的温度分布曲线图。

具体实施方式

本发明旨在提出一种产气井温度压力分布的数值模拟方法，对产气井的温度压力参数分 布情况进行准确模拟，从而为生产提供指导依据。

本发明中有可能涉及的参数释义如下：

A:内管面积，m²；

C_P:定压比热,J/Kg.K；

f:摩擦系数；

f(T_D):地层瞬态热传导时间方程；

g:重力加速度；

h:比焓；

k_ang:环中气体热传导性，J/m.K；

k_cem:座封热传导性，J/m.K；

k_e:地层热传导性，J/m.K；

l:井筒长度，l；

M:气体的平均摩尔重，g/mol；

P:流体压强，pa；

P_Pe:临界压强，pa；

rcem：座封外径，m

Re:雷诺数；

t:生产时间，s；

T:流体温度，K；

T_e:周围地层温度，K；

U_ti:总传热系数，W/m.K；

v:天然气速度，m/s；

w:总的质量流速；

z:压缩因子；

α：地层的热扩散率，m²/s。

下面结合附图及实施例对本发明的方案作更进一步的描述：

如图3所示，本发明产气井温度压力分布的数值模拟方法包括：

A、建立温度压力耦合微分方程组模型；

B、估计微分方程组相关参数；

C、使用差分法进行算法设计；

D、基于C#开发平台构建模拟软件；

E、在模拟软件中输入初始数据，获得温度压力分布模拟结果。

在具体实施上，步骤A中，基于假设的模型构建条件分别获得能量守恒方程、质量守恒 方程、动量守恒方程；

为符合实际产气井的情况，作如下假设：

(1)天然气的流动是单向稳定的；

(2)井筒的传热被认为是稳态的；

(3)基于无量纲瞬态导热时间方程，地层传热是不稳定的；

(4)导管和套管是同心的；

能量守恒方程：

$q+\frac{d}{dl}(h+v^2/2-gl\sin \theta )=0---\left(1\right)$

质量守恒方程：

$\rho \frac{dv}{dl}+v\frac{d\rho }{dl}=0---\left(2\right)$

动量守恒方程：

$-\frac{1}{\rho }\frac{dP}{dl}+f\frac{v^2}{2r_{ti}}+gsin\theta =\frac{vdv}{dl}---\left(3\right)$

耦合微分方程组模型：在微元dl上，从流体到座封-地层界面传导的辐射热为

$q=\frac{(-2\pi )r_{ti}{U}_{ti}}{w}(T-T_h)---\left(4\right)$

从座封-地层界面到周围地层传导的辐射热为

$q=\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\right](T-T_e)---\left(5\right)$

由(4)和(5)上，得到流体与周围地层热传导的微分方程：

$q=\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\right](T-T_e)---\left(6\right)$

将(1)和(6)整合，得到如下微分方程

$\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}\left(t_D\right)+ke}\right](T-T_e)+\frac{d}{dl}(h+v^2/2-gl\sin \theta )=0---\left(7\right)$

当流体在井筒中流动的时候，由于管径变化很小，Joule–Thomson系数可以忽略不计。因此， dh＝C_pdT，整体的焓变化量为

$\frac{2\pi }{w}\left[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\right](T-T_e)+C_p\frac{dT}{dl}+v\frac{dv}{dl}-g\sin \theta =0---\left(8\right)$

将(1)(2)(3)(8)放在一起，进行整理，可得到天然气井中的流体压力、温度、密度 和速率的微分组方程：

$\{\begin{array}{c}\frac{d\rho }{dl}=\frac{\frac{-8.314z\rho }{C_{p}M}[\frac{\left(2\pi \right)r_{ti}{U}_{ti}ke}{w[r_{ti}{U}_{ti}f\left(t_D\right)+ke]}(T-T_e)-g\cos \theta ]-f\frac{v^2\rho }{2r_{ti}}-g\rho \cos \theta }{v^2-(\frac{8.314{\mathrm{zv}}^2}{C_{p}M}+\frac{8.314Tz}{M})}\\ \frac{dv}{dl}=-\frac{v}{\rho }\frac{d\rho }{dz}\\ \frac{dP}{dl}=\rho g\cos \theta +f\frac{v^2\rho }{2r_{ti}}+v^2\frac{d\rho }{dl}\\ \frac{dT}{dl}=\frac{[\frac{v^2}{\rho }\frac{d\rho }{dl}+g\cos \theta -\frac{\left(2\pi \right)r_{ti}{U}_{ti}ke}{w[r_{ti}{U}_{ti}f\left(t_D\right)+ke]}(T-Te)]}{C_p}\end{array}---\left(9\right)$

用y_i(i＝1,2,3,4)代替ρ，v，P，T，则方程组可简化为

$\frac{{\mathrm{dy}}_i}{dl}=F_i(y_1,y_2,y_3,y_4),(i=1,2,3,4)---\left(10\right)$

步骤B中所述估计微分方程组相关参数包括：

B1.热传导系数估计：

$U_{ti}^{-1}=\frac{1}{hc+hr}+\frac{r_{ti}\ln \left(\frac{r_{cem}}{r_{co}}\right)}{k_{cem}}+\frac{r_{ti}\ln \left(\frac{r_{ci}}{r_{to}}\right)}{k_{ang}}---\left(11\right)$

上式表明在井筒系统中，套管/导管环隙间材料(气、油、水或者混合物)在决定热传导的 过程中起到重要的作用。

B2.计算定压比热C_P：

C_p＝1243+3.14T+7.931×10^-4T²-6.881×10^-7T³(12)

B3.地层瞬态热传导方程：

$\left(\begin{array}{cc}f\left(t_D\right)=1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}& (t_D\le 1.5)\\ f\left(t_D\right)=(1+\frac{0.6}{t_D})[0.4063+0.5\ln \left(t_D\right)]& (t_D>1.5)\end{array}\right)---\left(13\right)$

$t_D=\frac{t\alpha }{r_{wb}^2}$

B4.计算摩擦系数f：

$\frac{1}{\sqrt{f}}=1.14-2\lg (\frac{0.00001524}{r_{ti}}+\frac{21.25}{{\mathrm{Re}}^{0.9}})---\left(14\right).$

步骤C中使用差分法进行算法设计；由于C_P和θ不是常数，井筒必须划分为若干长度区间， 步长为h。假设井可以划分为n个区间，j(j＝1,2,…,n)为节点，n为井底。根据井底的温 度与压力，能够计算相应天然气的密度和流体速度。上述微分方程的边界条件为：

$P\left[n\right]=P_0,T\left[n\right]=T_0,\rho \left[n\right]=0.000001\times 3484.48{\gamma }_g\frac{P_0}{{\mathrm{zT}}_0},v\left[n\right]=\frac{101000\times 300000T_0}{293\times 86400P_{0}A}$

为了计算各节点的相关量，使用如下的迭代方程；

$y_i[j-1]=y_i\left[j\right]-\frac{h}{6}\left(a\right[i]+2b[i]+2c[i]+d[i\left]\right),(i=1,2,3,4;j=1,2,3,4...n)$

为使用以上方程，需要估计a[i],b[i],c[i]和d[i]:

$\left(\begin{array}{c}a\left[i\right]=F_i\left(y_1\right[n],y_2[n],y_3[n],y_4[n\left]\right)\\ b\left[i\right]=F_i\left(y_1\right[n]+\frac{h}{2}a[1],y_2[n]+\frac{h}{2}a[2],y_3[n]+\frac{h}{2}a[3],y_4[n]+\frac{h}{2}a[4\left]\right)\\ c\left[i\right]=F_i\left(y_1\right[n]+\frac{h}{2}b[1],y_2[n]+\frac{h}{2}b[2],y_3[n]+\frac{h}{2}b[3],y_4[n]+\frac{h}{2}b[4\left]\right)\\ d\left[i\right]=F_i\left(y_1\right[n]+hc[1],y_2[n]+hc[2],y_3[n]+hc[3],y_4[n]+hc[4\left]\right)\end{array}\right)$

步骤D中，基于C#开发平台构建模拟软件；软件结构如图4所示，其包括用户界面模块、 输入输出模块、算法模块和图标模块；

(1)用户界面模块：主要提供气井的初始信息，这些信息可划分为静态数据、流体数据和 生产数据：

(i)静态数据：描述静态下气井特征，简写为S。

(ii)流体数据：描述流体的状态，如水、油或气，简写为F。

(iii)生产数据：生产阶段给定的参数，简写为P。

(2)输入输出模块：该模块从存储器里读入、保存数据。使用MicrosoftSQLServer2000 建立井筒的数据库。

(3)算法模块：该模块将应用于所有的数学计算。能够计算出从井底到井口任意节点的温 度和压力。

(4)图表模块：该模块将数据图表化。通过CrystalReports可以画出温度压力分布的曲 线并输出结果。

在模拟软件设计完成后，就可以利用在模拟软件中输入初始数据，获得温度压力分布模 拟结果，具体算法流程如图5所示。

实施例：

对于某一口产气井，其输入数据见表1：

表1:产气井相关输入数据表

通过软件模拟，得到一系列结果。如图6和图7，温度和压力从井底到井口，越来越小。这 是因为重力、摩擦以及井筒对地层的热传导。温度在0到3000米，下降很多，但随后趋于稳 定，温度则随深度线性变化。图6显示，在同一井深，产气量越大，由于摩擦阻力额增加压 力越小，由于流速增大，温度升高。