CN103226641A - 深水气液两相流循环温度压力耦合计算方法 - Google Patents

Abstract

本发明公开了一种深水气液两相流循环温度压力耦合计算方法，包括以下步骤：1)计算网格垂向坐标；2)应用初始条件；3)从初始值出发，按照先钻柱内钻井液后环空钻井液的顺序迭代计算钻柱内和环空钻井液节点温度和压力数据，直至温度和压力都达到收敛条件，迭代结束，最后一次迭代计算结果为最终深海气液两相流井筒温度和压力模拟结果，保存并输出，该迭代命名为全局迭代。本发明的计算方法提高了计算精度。

Description

深水气液两相流循环温度压力耦合计算方法

技术领域

本发明涉及一种温度压力计算方法，具体涉及一种深水气液两相流循环温度压力耦合计算方法。

背景技术

随着欠平衡技术在陆地油田应用日趋成熟，近年来为了提高海上油井产能、降低钻井过程中对油气层的污染或解决井漏问题，已开始应用欠平衡钻井技术开发海洋油田。使用充气钻井液时，气液两相流钻井液形成的井筒温度场和压力场相互影响，属于耦合求解问题，井筒内气液两相流钻井液的温度和压力场数据是安全实施深水欠平衡钻井技术的关键。

目前，国内外学者对井筒传热问题的研究方法可以归纳为两大类：半瞬态法和全瞬态法。

半瞬态法认为井筒内钻井液传热速率远大于地层内岩石的传热速率，将井筒内传热视为稳态传热过程，地层内看作瞬态传热过程。该类方法的基础为Ramey模型，模型中将井筒看作是插入地层中的无限长线热源，依据井筒与地层的换热时间推导出了井筒与地层的换热量计算公式，适用于计算流体循环超过7天的井筒与地层换热问题。针对该模型缺陷，Jacques提出通过改进无因次函数f(t)的计算方法，使半稳态方法能适用于预测早期井筒与地层换热问题。

全瞬态法将井筒内换热和地层中换热均看作瞬态过程，原始模型由Raymond提出，其后，Keller在换热模型中加入了钻井液摩阻和机械能损失引起的内热源，David改进了离散方程组数值求解算法，加快了求解速度。

研究气液两相流循环压降以及界面含气率的模型可分为三类：均相流模型、分离流模型和基于流型的机理模型。国内外学者普遍认为分流型的机理模型能够更加准确的描述气液两相流在倾斜圆管中的真实流动状态，更适合计算井筒气液两相流压力降。

深水欠平衡钻井作业期间，由于气相的可压缩性，使得钻井液温度和压力与钻井液密度、流变性和热物性相互影响，海水区与地层区相互影响，钻柱内和环空温压场相互影响。而现有方法将气液两相流循环温度和压力分开单独求解，计算温度场时不考虑压力，计算压力时认为钻井液恒温，与深海钻井工况不符，不能用于模拟深水气液两相流温度压力场。

发明内容

为了解决上述技术问题，本发明提供一种深水气液两相流循环温度压力耦合计算方法。基于深水气液两相流流动、换热和压力传递特征，采用交错网格的全隐式有限体积法离散格式，耦合多换热区域，考虑钻井热源，温压对两相流钻井液热物性的影响，开发了本计算方法，并结合现场数据验证了本方法的有效性。

其特征在于，包括以下步骤：

1)节点划分

根据深水井筒井身结构和钻具结构，采用交错网格布置压力和温度节点，将温度节点布于网格控制体中心，压力节点布于网格控制体界面处。对每一个轴向网格都要分析和记录其所在轴向位置处的径向换热对象的几何和热物性信息。

2)确定网格垂向坐标

根据井眼轨迹确定网格中心的垂直坐标和网格垂直长度。

3)应用初始条件

海水区(泥线以上)钻柱内和环空所有节点的初始温度为节点垂直坐标对应深度处海水温度；地层区(泥线以下)钻柱内和环空所有节点的初始温度为网格中心垂直坐标对应垂直深度处地层原始温度；节点初始压力为气液两相钻井液静止状态下对应深度处的静液柱压力。

4)自上而下计算钻柱内气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下：

①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

②假定该节点的温度和压力等于上部相邻节点的温度和压力；

③取该节点与上部相邻节点压力的平均值为网格单元平均压力；

④计算环空气液两相流钻井液与钻柱外壁的强迫对流换热系数；

⑤计算气相和液相在网格中心温度和平均压力下的热物性参数；

⑥计算钻柱内气液两相流钻井液的压降梯度；

⑦计算钻柱内气液两相流与钻柱内壁的强迫对流换热系数；

⑧计算钻柱内钻井液与环空钻井液之间的热阻；

⑨计算钻柱内节点的新温度和压力；

⑩比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。

5)自下而上计算环空气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下：

①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

②假定该节点的温度和压力等于下部相邻节点的温度和压力；

③取该节点与下部相邻节点压力的平均值为网格单元平均压力；

④计算气相和液相在网格中心温度和平均压力下的热物性参数；

⑤计算环空气液两相流钻井液的压降梯度；

⑥计算环空气液两相流与钻柱外壁和环空壁面的强迫对流换热系数；

⑦计算环空钻井液与地层之间的热阻；

⑧计算环空节点的新温度和压力；

⑨比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。

6)比较节点新温度压力数据与Told和POld的差，确定相邻两次迭代所有节点的最大温度差值TDiffMax和压力差值PDiffMax，若最大温差和压力差满足收敛条件则终止迭代计算，保存计算结果，否则，转到步骤4重复计算直至获得收敛解。

一种深水气液两相流循环温度压力耦合计算方法的模拟器，由以下函数组成：

DataInput函数是TPWTP程序的数据输入函数，完成模拟井所有数据的输入，主要包括井身结构，钻具结构，海水深度，海水垂直温度分布，地层垂直温度分布，钻井液入口温度、排量、泵压、干度，气体种类，转速、扭矩，液体钻井液、钢材、水泥、地层隔水管绝热层等的参考状态热力学参数，地面温度；

TPField函数是TPWTP程序的总功能模块，按照程序流程图组装其它函数，完成深海气液两相流井筒温度和压力场计算及数据存储功能，计算钻柱内节点温度时需要调用TInDrillStem函数，计算环空节点温度则需要调用TAnnulus函数；

GridGeneration函数根据模拟井的井身结构、钻具结构、海水深度对温度和压力场求解域进行轴向分段和网格划分，存储网格节点的轴向几何信息及径向换热对象的几何和介质信息；

根据井眼轨迹计算网格中心的垂直坐标和网格垂直长度，函数功能要求调用DirectionParaCal函数根据井深和井斜角计算垂深；

TOriginGeneration函数根据海水和地层的垂直温度分布数据插值产生节点垂直深度处的原始温度；

Ini函数对求解域内网格节点变量应用初始条件，赋初值；

HTPPipe函数的功能是计算气液两相流管流强迫对流换热系数，为使程序可用于模拟气体钻井和液体钻井工况，HTPPipe还可通过调用HGasPipe和HLiquidPipe函数计算单气相或单液相的管流强迫对流换热系数；

HTPAnnulus函数的功能是计算环空气液两相流与内外壁的两个强迫对流换热系数，为使程序可用于模拟气体钻井和液体钻井工况，HTPAnnulus还可通过调用HGasAnnulus和HLiquidAnnulus函数计算单气相或单液相的环空强迫对流换热系数；

HSeaAcross函数用于计算海水横掠隔水管的强迫对流换热系数；

根据介质类型计算介质在给定温度和压力下的热力学性质的总集成函数，需要根据具体的介质类型空气、氮气、水、天然气、钻井液调用相应介质的热物性计算函数完成其功能；

DPDZ_BB：根据Beggs-Brill方法计算气液两相流压降梯度，实现该函数功能需要调用TPFriction_BB函数计算气液两相流摩阻压降；

DPDZ_HK：根据Hasan-Kabir方法计算气液两相流压降梯度，实现该函数功能需要调用TPFriction_HK函数计算气液两相流摩阻压降，

HeatResistance1：计算钻柱内钻井液到环空钻井液的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；

HeatResistance2：计算环空钻井液到地层的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；

TInDrillStem函数根据钻柱内气液两相钻井液能量守恒方程迭代计算钻柱内气液两相钻井液节点温度；

Tannulus：根据环空气液两相钻井液能量守恒方程迭代计算环空气液两相钻井液节点温度；

ThermalPhysics：计算气相在给定温度和压力下的密度和比焓，该函数功能需要根据气体类型调用空气；

EOS：根据气体类型调用相应的状态方程，根据温度和比体积计算气相的压力和比焓；

RKS：根据RKS模型计算给定温度和比体积时氮气的压力和比焓；

AirTP：计算空气在给定温度和比体积时氮气的压力和比焓，计算比焓时需要调用AirAOT、AirArT、AirArDen函数计算空气状态方程的相关偏导数；

AirAOT：计算空气的理想状态Helmholtz能对对比温度倒数的导数；

AirArT：计算空气余能对对比温度倒数的导数；

AirArDen：计算空气余能对对比密度的导数；

AirConductivity：计算空气给定温度和压力下的热导率，实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；

AirViscosity：计算空气给定温度和压力下的动力粘度，实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；

CpAir：计算空气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算空气的高温高压比焓；

CpNitrogen：计算氮气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算氮气的高温高压比焓，

NitrogenConductivity：计算氮气给定温度和压力下的热导率，实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度；

NitrogenViscosity：计算氮气给定温度和压力下的动力粘度，实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度，

本发明的有益效果：

(1)本发明的技术方案对钻柱内和环空气液两相钻井液按照有限体积法，对温度和压力耦合，海水区和地层区耦合，钻柱内和环空区域耦合，确定了算法的求解过程；

(2)本发明所述的计算方法是基于离散格式表示的深水气液两相流钻井液温度和压力耦合模型，采用交错网格的隐式有限差分法进行求解，提高了计算的稳定性；

(3)本发明考虑了温度和压力与钻井液密度及热传输特性的相互影响，钻头破岩生热，钻柱摩擦生热和钻井液摩阻生热等的热源作用，提高了计算的准确性；

(4)利用本发明技术所述方法的计算数据与现场实例数据进行了对比，验证了该模拟器的准确性，计算误差不超过5％，可以投入现场使用。

附图说明

图1为半潜式钻井平台作业时的物理模型图；

图2为钻柱内和环空温度及压力节点划分图；

图3为本发明方法计算总流程图；

图4为本发明方法钻柱内温度和压力计算流程图；

图5为本发明方法环空内温度和压力计算流程图；

图6为本发明TPWTP程序主要模块结构图。

具体实施方式

下面结合附图具体实施方式对本发明的方法作进一步详细地说明。

深水钻井作业时的物理模型如图1所示，上部被海水包围，下部被地层包围。气液两相流在井筒内循环时，依据传热过程，将模型求解区域划分为5个区：①钻柱内流体区，两相流钻井液自钻井泵流入钻柱内，沿钻柱向下一直到井底；②钻柱与井壁之间的环空区，两相流钻井液从井底进入环空，向上流动，直至井口；③井筒/地层(海水)界面区，井筒与地层和海水的交界面；④地层区；⑤海水区。

深水气液两相流循环温度压力场计算步骤

1)按交错网格布置节点

在图1所示的深海井筒传热物理模型所划分的5个计算区域内，将节点布置于钻柱内和环空两个区域内。考虑到压力求解的稳定和收敛问题，采用交错网格布置压力和温度节点，将温度节点布于网格控制体中心，压力节点布于网格控制体界面处。对每一个轴向网格都要分析和记录其所在轴向位置处的径向换热对象的几何和热物性信息。节点轴向序号自井口向井底递增，井口节点轴向序号为0。

(1)轴向网格划分

轴向网格划分时，将求解区域依据隔水管深度、井身结构和钻具结构进行轴向分段，分段原则为：每一段内涉及到的所有传热对象的几何尺寸只有一种，分段信息记录在Sections数组中。然后自上至下对每一段再根据段长进行轴向网格划分，每一段内的网格控制体积相等。算法中使用NodeZ和NodeDZ两个数组记录网格中心轴向坐标和网格轴向长度。为节省程序内存，使用NodeInSection数组记录网格所在的段索引，大大减小了节点径向信息存储需要开支的存储空间。

(2)径向信息记录

从上到下通过扫描钻具结构和井身结构数据，分析每一段涉及到的径向传热对象几何信息和介质种类，记录在Annulus二维数组中。假如某段有两层套管，套管外均为水泥环，则井筒内径向传热对象有：最内层为钻柱内钻井液，往外依次为钻柱、钻柱外环空钻井液、最内层套管I、套管I外水泥环、外层套管II、套管II外水泥环，需要记录的信息有：钻柱内钻井液的介质类型(钻井液)、钻柱内钻井液的径向尺寸、钻柱管体介质类型(钢材)、钻柱管体径向尺寸、环空钻井液介质类型(钻井液)、环空钻井液径向尺寸、套管I介质(钢材)、套管I径向尺寸、套管I环空介质类型(水泥石)、套管I环空径向几何尺寸、套管II介质类型(钢材)、套管II径向几何尺寸、套管II环空介质类型(水泥石)、套管II环空径向几何尺寸。

2)确定网格垂向坐标

根据井眼轨迹确定网格中心的垂直坐标和网格垂直长度。

3)应用初始条件

海水区(泥线以上)钻柱内和环空所有节点的初始温度为节点垂直坐标对应深度处海水温度；地层区(泥线以下)钻柱内和环空所有节点的初始温度为网格中心垂直坐标对应垂直深度处地层原始温度；节点初始压力为气液两相钻井液静止状态下对应深度处的静液柱压力。

4)自上而下计算钻柱内气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下：

①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

②假定该节点的温度和压力等于上部相邻节点的温度和压力；

③取该节点与上部相邻节点压力的平均值为网格单元平均压力；

④确定环空气液两相流钻井液与钻柱外壁的强迫对流换热系数w₂；

⑤确定气相和液相在网格中心温度和平均压力下的热物性参数：密度、定压比热、热导率、动力粘度等；

气相要根据具体气体类型是空气还是氮气，选择合适的高温高压热物性参数，液相则需根据钻井液类型选择合适的热物性计算方法。

⑥按式(1)确定钻柱内气液两相流钻井液的压降梯度(dp/dZ)_i；

${\left(\frac{\mathrm{dp}}{\mathrm{dZ}}\right)}_i=-\frac{\∂\left({\mathrm{\ρ}}_1{v}_1^2\right)}{\∂Z}+{\mathrm{\ρ}}_{1}g\cos \mathrm{\θ}-\frac{{2\mathrm{\ρ}}_1{f}_1{v}_1^2}{d_1}---\left(1\right)$

式中f₁为钻柱内气液两相流摩阻系数，基于两相流力学模型判别方法，根据节点处的流动状态(泡状流、段塞流、块状流和环状流)和截面含汽率来确定。

⑦确定钻柱内气液两相流与钻柱内壁的强迫对流换热系数w₁；

⑧按式(2)计算钻柱内钻井液与环空钻井液之间的热阻R₁₂；

钻柱内流体与环空流体换热过程中发生的热交换有：钻柱内流体与钻柱内壁的对流换热、钻柱管体的径向热传导、钻柱外壁与环空流体的对流换热。根据热通量相等，得总热阻表达式：

$R_{12}=\frac{1}{{\mathrm{\πd}}_1{w}_1}+\frac{1}{{\mathrm{\πd}}_2{w}_2}+\frac{\ln (d_2/d_1)}{{2\mathrm{\πk}}_{\mathrm{dp}}}---\left(2\right)$

式中：d₁为钻柱内径，m；d₂为钻柱外径，m；w₁为钻柱内钻井液与钻柱内壁的对流换热系数，W/(m2·K)；w₂为环空钻井液与钻柱外壁的对流换热系数，W/(m2·K)；k_dp为钻柱热导率，W/(m·K)。

⑨按式(3)计算钻柱内节点的新温度，按式(5)计算节点新压力；

钻柱内节点i的温度控制离散方程： $\frac{(T_{2,i}-T_{1,i})}{R_{12,i}}{\mathrm{\Δz}}_i{-\mathrm{mc}}_p(T_{1,i}-T_{1,i-1})+\mathrm{mg}\cos \mathrm{\θ\Δ}{z}_i$

$-m(\frac{v_{1,i+1/2}^2}{2}-\frac{v_{1,i-1/2}^2}{2})+{\mathrm{mh}}_p\left(p_{1,i+1/2}{-p}_{1,i-1/2}\right)+{\stackrel{\·}{q}}_{1,i}{\mathrm{\Δz}}_i=0---\left(3\right)$

钻柱内节点热源由气液两相钻井液流动摩阻产生，

${\stackrel{\·}{q}}_f=Q\×{\left(\frac{\mathrm{dp}}{\mathrm{dZ}}\right)}_f---\left(4\right)$

式中：T为温度，℃；p为压力，Pa；k为热导率，W/(m·K)；ρ为气液两相根据气相体积分数折算出的混合密度，kg/m³；v为气液两相根据气相体积分数折算出的钻井液流速，m/s；

为单位长度生热量，W/m；Z为深度坐标，m；r为径向坐标，m；R为总热阻，K/W；

为气液两相钻井液定温条件下比焓随压力的变化；m为气液两相钻井液质量流量，kg/s；g为重力加速度，m/s²；θ为井斜角，℃。下标：1为钻柱内钻井液；2为环空；3在地层段表示井筒与地层交界面，在海水段表示隔水管与海水交界面；e为地层。

钻柱内节点i压力： $p_i=p_{i-1}+{\left(\frac{\mathrm{dp}}{\mathrm{dZ}}\right)}_i---\left(5\right)$

⑩比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。

5)自下而上计算环空气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下：

①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

②假定该节点的温度和压力等于下部相邻节点的温度和压力；

③取该节点与下部相邻节点压力的平均值为网格单元平均压力；

④计算气相和液相在网格中心温度和平均压力下的热物性参数；

⑤计算环空气液两相流钻井液的压降梯度(dp/dZ)_i；

${\left(\frac{\mathrm{dp}}{\mathrm{dZ}}\right)}_i=-\frac{\∂\left({\mathrm{\ρ}}_3{v}_3^2\right)}{\∂Z}+{\mathrm{\ρ}}_{3}g\cos \mathrm{\θ}+\frac{{2\mathrm{\ρ}}_3{f}_3{v}_3^2}{(d_3-d_2)}---\left(6\right)$

⑥计算环空气液两相流与钻柱外壁和环空壁面的强迫对流换热系数w₂和w₃；

⑦计算环空钻井液与环境(地层或海水)之间的热阻R₃；

a对地层段节点

$R_3=R_{3w}+\frac{f\left(t\right)}{{2\mathrm{\πk}}_e}---\left(7\right)$

R_3w为环空流体与井筒/地层界面换热的总热阻。环空流体与井筒/地层界面换热过程中发生的热交换有：环空内流体与井筒壁面的对流换热、套管管体的径向热传导、相邻套管柱形成的环空间隙内的水泥石的径向热传导。

$R_{3w}=\frac{1}{{2\mathrm{\πr}}_3{w}_3}+\underset{j=1}{\overset{N_c}{\mathrm{\Σ}}}\frac{\ln (r_{\mathrm{jO}}/r_{\mathrm{jI}})}{{2\mathrm{\πk}}_c}+\underset{s=1}{\overset{N_c}{\mathrm{\Σ}}}\frac{\ln (r_{\mathrm{sO}}/r_{\mathrm{sI}})}{{2\mathrm{\πk}}_s}---\left(8\right)$

式中：Nc为环空流体外环形介质层数，包括套管、套管间水泥环或钻井液；下标c表示水泥环，s表示钢材。

b对海水段节点

$R_3=R_{3r}+\frac{1}{{2\mathrm{\πr}}_{\mathrm{ro}}{w}_r}---\left(9\right)$

式中R_3r为环空流体与井筒/海水界面换热的总热阻。环空流体与井筒/海水界面换热过程中发生的热交换有：环空内流体与隔水管壁面的对流换热、隔水管管体的径向热传导。

$R_{3r}=\frac{1}{{\mathrm{\πd}}_{r1}{w}_3}+\frac{\ln (d_{r0}/d_{r1})}{{2\mathrm{\πk}}_r}---\left(10\right)$

式中：d_ri为隔水管内径，m；d_ro为隔水管外径，m；k_r为隔水管热导率，W/(m·K)。

⑧计算环空节点的新温度(式11)和压力(式15)；

环空节点i的温度控制离散方程：

对地层段环空节点有，

$(T_{1,i}-T_{3,i}){\mathrm{\Δz}}_i/R_{13,i}+{\mathrm{mc}}_p(T_{3,i+1}-T_{3,i})+(T_{e,i}-T_{3,i}){\mathrm{\Δz}}_i/R_{3e,i}-\mathrm{mg}\cos {\mathrm{\θ\Δz}}_i$

$+m(\frac{v_{3,i+1/2}^2}{2}-\frac{v_{3,i-1/2}^2}{2})+{\mathrm{mh}}_p(p_{3,i+1/2}-p_{3,i-1/2})+{\stackrel{\·}{q}}_{3,i}{\mathrm{\Δz}}_i=0---\left(11\right)$

对海水段环空节点有，

$(T_{1,i}-T_{3,i}){\mathrm{\Δz}}_i/R_{13,i}+{\mathrm{mc}}_p(T_{3,i+1}-T_{3,i})+(T_{w,i}-T_{3,i})\mathrm{\Δ}/R_{3w,i}-\mathrm{mg}\cos {\mathrm{\θ\Δz}}_i$

$+m(\frac{v_{3,i+1/2}^2}{2}-\frac{v_{3,i-1/2}^2}{2})+{\mathrm{mh}}_p(p_{3,i+1/2}-p_{3,i-1/2})+{\stackrel{\·}{q}}_{3,i}{\mathrm{\Δz}}_i=0---\left(12\right)$

环空内节点热源确定。

环空内节点热源(非钻头处节点)由两部分组成：气液两相钻井液流动摩阻生热

(按式4计算)和钻柱与井壁摩擦生热

${\stackrel{\·}{q}}_s=\mathrm{\β}\×Z_j\×2\mathrm{\πM}\×\frac{\mathrm{RPM}}{60D}---\left(13\right)$

对钻头处节点，热源由钻头破岩摩擦生热产生，

${\stackrel{\·}{q}}_b=\mathrm{\α}\×2\mathrm{\πM}\×\frac{\mathrm{RPM}}{60}---\left(14\right)$

环空内节点i压力： $p_i=p_{i+1}-{\left(\frac{\mathrm{dp}}{\mathrm{dZ}}\right)}_i---\left(15\right)$

⑨比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。

6)比较节点新温度压力数据与Told和POld的差，确定相邻两次迭代所有节点的最大温度差值TDiffMax和压力差值PDiffMax，若最大温差和压力差满足收敛条件则终止迭代计算，保存计算结果，否则，转到步骤4)重复计算直至获得收敛解。

有效性分析

根据专利所述方法与Keller井筒温度场数据和Louisiana州立大学的全尺寸两相流压力实验数据进行了对比。

与Keller温度数据对比

利用Keller井基本数据，计算了钻井液循环6、12、18、24、144小时后的井底温度，与Keller解对比见表1。

表1与Keller井井底温度对比

Table1 Comparisons of predicted bottom-hole temperatures on the Keller well

表1中数据表明本方法与Keller的温度数据吻合程度较高，最大相对误差位2.59％。

气液两相流钻井液循环压力验证

根据Louisiana州立大学1#钻井液的注氮气两相流循环压力实验，利用本文模型和程序对两相流循环压力进行了预测。实验条件：管路内径：50.673mm；管路长度：902.82m；管路状态：垂直管；流动方向：垂直向上流；钻井液密度：1030.50kg/m³；钻井液塑性粘度：0.004Pa.s；钻井液动切力：0.48Pa；流性指数：0.65；稠度系数：0.061Pa.s；注气量：306.76l/s；室温：23.89℃；地温梯度2.00℃/100m。本文预测值与实验数据对比见表2。

表2气液两相流钻井液循环压力计算结果与实验结果对比

Table2 Comparisons of gas-liquid flow pressure drop between predicted andexperimental values

表2数据表明，本方法计算的两相流压力与实验数据的相对误差小于5％，与现有技术的两相流压力计算模型精度相当，甚至略高，满足工程计算需要。

根据气液两相流深海井筒温度和压力分布计算方法开发了形成了TPWTP模拟器。

TPWTP程序编写方式同WHT程序，共由57个函数构成，其主要功能模块(函数)如下，具体参见图6.

(1)DataInput

DataInput函数是TPWTP程序的数据输入函数，完成模拟井所有数据的输入，主要包括井身结构，钻具结构，海水深度，海水垂直温度分布，地层垂直温度分布，钻井液入口温度、排量、泵压、干度，气体种类，转速、扭矩，液体钻井液、钢材、水泥、地层隔水管绝热层等的参考状态热力学参数，地面温度等。

(2)TPField

TPField函数是TPWTP程序的总功能模块，按照程序流程图组装其它函数，完成深海气液两相流井筒温度和压力场计算及数据存储功能。计算钻柱内节点温度时需要调用TInDrillStem函数，计算环空节点温度则需要调用TAnnulus函数。

(3)GridGeneration

GridGeneration函数根据模拟井的井身结构、钻具结构、海水深度对温度和压力场求解域进行轴向分段和网格划分，存储网格节点的轴向几何信息及径向换热对象的几何和介质信息(网格控制体内的介质有：气液两相钻井液、钢材、水泥、地层、隔水管绝热层等五种介质类型，每种介质都有相应的热力学物理性质)。

(4)VerticalCoordinate

根据井眼轨迹计算网格中心的垂直坐标和网格垂直长度。函数功能要求调用DirectionParaCal函数根据井深和井斜角计算垂深。

(5)TOriginGeneration

TOriginGeneration函数根据海水和地层的垂直温度分布数据插值产生节点垂直深度处的原始温度。

(6)Ini

Ini函数对求解域内网格节点变量应用初始条件，赋初值。

(7)HTPPipe

HTPPipe函数的功能是计算气液两相流管流强迫对流换热系数。为使程序可用于模拟气体钻井和液体钻井工况，HTPPipe还可通过调用HGasPipe和HLiquidPipe函数计算单气相或单液相的管流强迫对流换热系数。

(13)HTPAnnulus

HTPAnnulus函数的功能是计算环空气液两相流与内外壁的两个强迫对流换热系数。为使程序可用于模拟气体钻井和液体钻井工况，HTPAnnulus还可通过调用HGasAnnulus和HLiquidAnnulus函数计算单气相或单液相的环空强迫对流换热系数。

(14)HSeaAcross

HSeaAcross函数用于计算海水横掠隔水管的强迫对流换热系数。

(15)PPTP

根据介质类型计算介质在给定温度和压力下的热力学性质的总集成函数，需要根据具体的介质类型(空气、氮气、水、天然气、钻井液)调用相应介质的热物性计算函数完成其功能。

(16)DPDZ_BB

根据Beggs-Brill方法计算气液两相流压降梯度。实现该函数功能需要调用TPFriction_BB函数计算气液两相流摩阻压降。

(17)DPDZ_HK

根据Hasan-Kabir方法计算气液两相流压降梯度。实现该函数功能需要调用TPFriction_HK函数计算气液两相流摩阻压降。

(18)HeatResistance1

计算钻柱内钻井液到环空钻井液的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数。

(19)HeatResi stance2

计算环空钻井液到地层的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数。

(20)TInDrillStem

TInDrillStem函数根据钻柱内气液两相钻井液能量守恒方程迭代计算钻柱内气液两相钻井液节点温度。

(21)TAnnulus

根据环空气液两相钻井液能量守恒方程迭代计算环空气液两相钻井液节点温度。

(22)ThermalPhysics

计算气相在给定温度和压力下的密度和比焓。该函数功能需要根据气体类型调用空气或氮气的状态方程计算特定温度和比体积下的压力，然后利用数值方法计算气相的密度和比焓。

(23)EOS

根据气体类型调用相应的状态方程，根据温度和比体积计算气相的压力和比焓。

(24)RKS

根据RKS模型计算给定温度和比体积时氮气的压力和比焓。

(25)AirTP

计算空气在给定温度和比体积时氮气的压力和比焓，计算比焓时需要调用AirAOT、AirArT、AirArDen等函数计算空气状态方程的相关偏导数。

(26)AirAOT

计算空气的理想状态Helmholtz能对对比温度倒数的导数。

(27)AirArT

计算空气余能对对比温度倒数的导数。

(28)AirArDen

计算空气余能对对比密度的导数。

(29)AirConductivity

计算空气给定温度和压力下的热导率。实现函数功能需要调用ThermalPhysics函数计算空气的真实密度。

(30)AirViscosity

计算空气给定温度和压力下的动力粘度。实现函数功能需要调用ThermalPhysics函数计算空气的真实密度。

(31)CpAir

计算空气给定温度和压力下的定压比热。实现函数功能需要调用ThermalPhysics函数计算空气的高温高压比焓。

(32)CpNitrogen

计算氮气给定温度和压力下的定压比热。实现函数功能需要调用ThermalPhysics函数计算氮气的高温高压比焓。

(33)NitrogenConductivity

计算氮气给定温度和压力下的热导率。实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度。

(34)NitrogenViscosity

计算氮气给定温度和压力下的动力粘度。实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。

Claims (4)

1.一种深水气液两相流循环温度压力耦合计算方法，其特征在于，包括以下步骤：

1)计算网格垂向坐标

根据井眼轨迹计算网格中心的垂直坐标和网格垂直长度；

2)应用初始条件

海水区钻柱内和环空所有节点的初始温度为节点垂直坐标对应深度处海水温度，海水温度可由程序使用人员按水深输入，也可根据深海温度垂直分布模型考虑季节因素计算得到；

地层区钻柱内和环空所有节点的初始温度为网格中心垂直坐标对应垂直深度处地层原始温度，该原始温度同样可由用户输入也可根据地温梯度模型计算得到；

3)从初始值出发，按照先钻柱内钻井液后环空钻井液的顺序迭代计算钻柱内和环空钻井液节点温度和压力数据，直至温度和压力都达到收敛条件，迭代结束，最后一次迭代计算结果为最终深海气液两相流井筒温度和压力模拟结果，保存并输出，该迭代命名为全局迭代，由于对钻柱内钻井液和环空钻井液每个节点求解温度和压力值时都存在温度、压力与气液两相钻井液热物性的相互影响，故对每一个节点的求解也需要迭代计算，直至取得稳定收敛解，该迭代命名为子迭代。

2.根据权利要求1所述的深水气液两相流循环温度压力耦合计算方法，其特征在于，所述步骤3)中每一次全局迭代时钻柱内钻井液某节点的温度和压力计算步骤如下：

A用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

B假定该节点的温度和压力等于上部相邻节点的温度和压力；

C取该节点与上部相邻节点压力的平均值为网格单元平均压力；

D计算环空气液两相流钻井液与钻柱外壁的强迫对流换热系数；

E计算气相和液相在网格中心温度和平均压力下的热物性参数；

F计算钻柱内气液两相流钻井液的压降梯度；

G计算钻柱内气液两相流与钻柱内壁的强迫对流换热系数；

H计算钻柱内钻井液与环空钻井液之间的热阻；

I计算钻柱内节点的新温度和压力；

J比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤C，重复执行，直至获得收敛解，作为该次迭代的最终解。

3.根据权利要求1所述的深水气液两相流循环温度压力耦合计算方法，其特征在于，所述步骤3)中每一次全局迭代时环空钻井液某节点的温度和压力计算步骤如下：

K用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；

L假定该节点的温度和压力等于下部相邻节点的温度和压力；

M取该节点与下部相邻节点压力的平均值为网格单元平均压力；

N计算气相和液相在网格中心温度和平均压力下的热物性参数；

O计算环空气液两相流钻井液的压降梯度；

P计算环空气液两相流与钻柱外壁和环空壁面的强迫对流换热系数；

Q计算环空钻井液与地层之间的热阻；

R计算环空节点的新温度和压力；

S比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤M，重复执行，直至获得收敛解，作为该次迭代的最终解。

T比较节点新温度压力数据与Told和POld的差，确定相邻两次迭代所有节点的最大温度差值TDiffMax和压力差值PDiffMax，以判断全局迭代计算是否满足结束条件。

4.一种实现权利要求1所述的深水气液两相流循环温度压力耦合计算方法的模拟器，其特征在于，由以下函数组成：

DataInput函数是TPWTP程序的数据输入函数，完成模拟井所有数据的输入，主要包括井身结构，钻具结构，海水深度，海水垂直温度分布，地层垂直温度分布，钻井液入口温度、排量、泵压、干度，气体种类，转速、扭矩，液体钻井液、钢材、水泥、地层隔水管绝热层等的参考状态热力学参数，地面温度；

TPField函数是TPWTP程序的总功能模块，按照程序流程图组装其它函数，完成深海气液两相流井筒温度和压力场计算及数据存储功能，计算钻柱内节点温度时需要调用TInDrillStem函数，计算环空节点温度则需要调用TAnnulus函数；

GridGeneration函数根据模拟井的井身结构、钻具结构、海水深度对温度和压力场求解域进行轴向分段和网格划分，存储网格节点的轴向几何信息及径向换热对象的几何和介质信息；

根据井眼轨迹计算网格中心的垂直坐标和网格垂直长度，函数功能要求调用DirectionParaCal函数根据井深和井斜角计算垂深；

TOriginGeneration函数根据海水和地层的垂直温度分布数据插值产生节点垂直深度处的原始温度；

Ini函数对求解域内网格节点变量应用初始条件，赋初值；

HTPPipe函数的功能是计算气液两相流管流强迫对流换热系数，为使程序可用于模拟气体钻井和液体钻井工况，HTPPipe还可通过调用HGasPipe和HLiquidPipe函数计算单气相或单液相的管流强迫对流换热系数；

HTPAnnulus函数的功能是计算环空气液两相流与内外壁的两个强迫对流换热系数，为使程序可用于模拟气体钻井和液体钻井工况，HTPAnnulus还可通过调用HGasAnnulus和HLiquidAnnulus函数计算单气相或单液相的环空强迫对流换热系数；

HSeaAcross函数用于计算海水横掠隔水管的强迫对流换热系数；

根据介质类型计算介质在给定温度和压力下的热力学性质的总集成函数，需要根据具体的介质类型空气、氮气、水、天然气、钻井液调用相应介质的热物性计算函数完成其功能；

DPDZ_BB：根据Beggs-Brill方法计算气液两相流压降梯度，实现该函数功能需要调用TPFriction_BB函数计算气液两相流摩阻压降；

DPDZ_HK：根据Hasan-Kabir方法计算气液两相流压降梯度，实现该函数功能需要调用TPFriction_HK函数计算气液两相流摩阻压降，

HeatResistance1：计算钻柱内钻井液到环空钻井液的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；

HeatResistance2：计算环空钻井液到地层的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；

TInDrillStem函数根据钻柱内气液两相钻井液能量守恒方程迭代计算钻柱内气液两相钻井液节点温度；

Tannulus：根据环空气液两相钻井液能量守恒方程迭代计算环空气液两相钻井液节点温度；

ThermalPhysics：计算气相在给定温度和压力下的密度和比焓，该函数功能需要根据气体类型调用空气；

EOS：根据气体类型调用相应的状态方程，根据温度和比体积计算气相的压力和比焓；

RKS：根据RKS模型计算给定温度和比体积时氮气的压力和比焓；

AirTP：计算空气在给定温度和比体积时氮气的压力和比焓，计算比焓时需要调用AirAOT、AirArT、AirArDen函数计算空气状态方程的相关偏导数；

AirAOT：计算空气的理想状态Helmholtz能对对比温度倒数的导数；

AirArT：计算空气余能对对比温度倒数的导数；

AirArDen：计算空气余能对对比密度的导数；

AirConductivity：计算空气给定温度和压力下的热导率，实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；

AirViscosity：计算空气给定温度和压力下的动力粘度，实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；

CpAir：计算空气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算空气的高温高压比焓；

CpNitrogen：计算氮气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算氮气的高温高压比焓，

NitrogenConductivity：计算氮气给定温度和压力下的热导率，实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度；

NitrogenViscosity：计算氮气给定温度和压力下的动力粘度，实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度。

Images