专利名称：深水气液两相流循环温度压力耦合计算方法
技术领域：
本发明涉及一种温度压力计算方法，具体涉及一种深水气液两相流循环温度压力耦合计算方法。
背景技术：
随着欠平衡技术在陆地油田应用日趋成熟，近年来为了提高海上油井产能、降低钻井过程中对油气层的污染或解决井漏问题，已开始应用欠平衡钻井技术开发海洋油田。使用充气钻井液时，气液两相流钻井液形成的井筒温度场和压力场相互影响，属于耦合求解问题，井筒内气液两相流钻井液的温度和压力场数据是安全实施深水欠平衡钻井技术的关键。目前，国内外学者对井筒传热问题的研究方法可以归纳为两大类:半瞬态法和全瞬态法。半瞬态法认为井筒内钻井液传热速率远大于地层内岩石的传热速率，将井筒内传热视为稳态传热过程，地层内看作瞬态传热过程。该类方法的基础为Ramey模型，模型中将井筒看作是插入地层中的无限长线热源，依据井筒与地层的换热时间推导出了井筒与地层的换热量计算公式，适用于计算流体循环超过7天的井筒与地层换热问题。针对该模型缺陷，Jacques提出通过改进无因次函数f(t)的计算方法，使半稳态方法能适用于预测早期井筒与地层换热问题。全瞬态法将井筒内换热 和地层中换热均看作瞬态过程，原始模型由Raymond提出，其后，KelIer在换热模型中加入了钻井液摩阻和机械能损失弓I起的内热源，David改进了离散方程组数值求解算法，加快了求解速度。研究气液两相流循环压降以及界面含气率的模型可分为三类:均相流模型、分离流模型和基于流型的机理模型。国内外学者普遍认为分流型的机理模型能够更加准确的描述气液两相流在倾斜圆管中的真实流动状态，更适合计算井筒气液两相流压力降。深水欠平衡钻井作业期间，由于气相的可压缩性，使得钻井液温度和压力与钻井液密度、流变性和热物性相互影响，海水区与地层区相互影响，钻柱内和环空温压场相互影响。而现有方法将气液两相流循环温度和压力分开单独求解，计算温度场时不考虑压力，计算压力时认为钻井液恒温，与深海钻井工况不符，不能用于模拟深水气液两相流温度压力场。
发明内容
为了解决上述技术问题，本发明提供一种深水气液两相流循环温度压力耦合计算方法。基于深水气液两相流流动、换热和压力传递特征，采用交错网格的全隐式有限体积法离散格式，耦合多换热区域，考虑钻井热源，温压对两相流钻井液热物性的影响，开发了本计算方法，并结合现场数据验证了本方法的有效性。其特征在于，包括以下步骤:
I)节点划分根据深水井筒井身结构和钻具结构，采用交错网格布置压力和温度节点，将温度节点布于网格控制体中心，压力节点布于网格控制体界面处。对每一个轴向网格都要分析和记录其所在轴向位置处的径向换热对象的几何和热物性信息。2)确定网格垂向坐标根据井眼轨迹确定网格中心的垂直坐标和网格垂直长度。3)应用初始条件海水区(泥线以上)钻柱内和环空所有节点的初始温度为节点垂直坐标对应深度处海水温度；地层区(泥线以下)钻柱内和环空所有节点的初始温度为网格中心垂直坐标对应垂直深度处地层原始温度；节点初始压力为气液两相钻井液静止状态下对应深度处的静液柱压力。4)自上而下计算钻柱内气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下:①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；②假定该节点的温度和压力等于上部相邻节点的温度和压力；③取该节点与上部相邻节点压力的平均值为网格单元平均压力；④计算环空气液两相流钻井液与钻柱外壁的强迫对流换热系数；⑤计算气相和液相在网格中心温度和平均压力下的热物性参数；
⑥计算钻柱内气液两相流钻井液的压降梯度；⑦计算钻柱内气液两相流与钻柱内壁的强迫对流换热系数；⑧计算钻柱内钻井液与环空钻井液之间的热阻；⑨计算钻柱内节点的新温度和压力；⑩比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。5)自下而上计算环空气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下:①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；②假定该节点的温度和压力等于下部相邻节点的温度和压力；③取该节点与下部相邻节点压力的平均值为网格单元平均压力；④计算气相和液相在网格中心温度和平均压力下的热物性参数；⑤计算环空气液两相流钻井液的压降梯度；⑥计算环空气液两相流与钻柱外壁和环空壁面的强迫对流换热系数；⑦计算环空钻井液与地层之间的热阻；⑧计算环空节点的新温度和压力；⑨比较并记录本次迭代节点初始温度和压力与新计算出的节点温度和压力的差值，若达到收敛条件，则该节点本次温度和压力迭代计算结束，否则，以新计算出的温度和压力作为初始值，再转到步骤③，重复执行，直至获得收敛解，作为该次迭代的最终解。6)比较节点新温度压力数据与Told和POld的差，确定相邻两次迭代所有节点的最大温度差值TDiffMax和压力差值H)iffMax，若最大温差和压力差满足收敛条件则终止迭代计算，保存计算结果，否则，转到步骤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方法计算气液两相流压降梯度,实现该函数功能需要调用TPFrictionJlK函数计算气液两相流摩阻压降，HeatResistancel:计算钻柱内钻井液到环空钻井液的换热热阻,实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；HeatResistancd:计算环空钻井液到地层的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；TInDrillStem函数根据钻柱内气液两相钻井液能量守恒方程迭代计算钻柱内气液两相钻井液节点温度；Tannulus:根据环空气液两相钻井液能量守恒方程迭代计算环空气液两相钻井液T1点温度；ThermalPhysics:计算气相在给定温度和压力下的密度和比j:含,该函数功能需要根据气体类型调用空气；E0S:根据气体类型调用相应的状态方程，根据温度和比体积计算气相的压力和比焓；RKS:根据RKS模型计算给定温度和比体积时氮气的压力和比焓；AirTP:计算空气在给定温度和比体积时氮气的压力和比焓，计算比焓时需要调用AirAOT、AirArT、AirArDen函数计算空气状态方程的相关偏导数；AirAOT:计算空气的理想状态Helmholtz能对对比温度倒数的导数；AirArT:计算空气余能对对比温度倒数的导数；AirArDen:计算空气余能对对比密度的导数；AirConductivity:计算空气给定温度和压力下的热导率,实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；AirViscosity:计算空气给定温度和压力下的动力粘度,实现函数功能需要调用ThermalPhysics函数计算空气的真实密度；CpAir:计算空气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算空气 的高温高压比j:含；CpNitrogen:计算氮气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算氮气的高温高压比j:含，NitrogenConductivity:计算氮气给定温度和压力下的热导率,实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度；NitrogenViscosity:计算氮气给定温度和压力下的动力粘度,实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度，本发明的有益效果:(I)本发明的技术方案对钻柱内和环空气液两相钻井液按照有限体积法，对温度和压力耦合，海水区和地层区耦合，钻柱内和环空区域耦合，确定了算法的求解过程；(2)本发明所述的计算方法是基于离散格式表示的深水气液两相流钻井液温度和压力耦合模型，采用交错网格的隐式有限差分法进行求解，提高了计算的稳定性；(3)本发明考虑了温度和压力与钻井液密度及热传输特性的相互影响，钻头破岩生热，钻柱摩擦生热和钻井液摩阻生热等的热源作用，提高了计算的准确性；(4)利用本发明技术所述方法的计算数据与现场实例数据进行了对比，验证了该模拟器的准确性，计算误差不超过5 %，可以投入现场使用。
图1为半潜式钻井平台作业时的物理模型图；图2为钻柱内和环空温度及压力节点划分图；图3为本发明方法计算总流程图4为本发明方法钻柱内温度和压力计算流程图；图5为本发明方法环空内温度和压力计算流程图；图6为本发明TPWTP程序主要模块结构图。
具体实施例方式下面结合附图具体实施方式
对本发明的方法作进一步详细地说明。深水钻井作业时的物理模型如图1所示，上部被海水包围，下部被地层包围。气液两相流在井筒内循环时，依据传热过程，将模型求解区域划分为5个区:①钻柱内流体区，两相流钻井液自钻井泵流入钻柱内，沿钻柱向下一直到井底钻柱与井壁之间的环空区，两相流钻井液从井底进入环空，向上流动，直至井口 ；③井筒/地层(海水)界面区，井筒与地层和海水的交界面地层区；⑤海水区。深水气液两相流循环温度压力场计算步骤I)按交错网格布置节点
在图1所示的深海井筒传热物理模型所划分的5个计算区域内，将节点布置于钻柱内和环空两个区域内。考虑到压力求解的稳定和收敛问题，采用交错网格布置压力和温度节点，将温度节点布于网格控制体中心，压力节点布于网格控制体界面处。对每一个轴向网格都要分析和记录其所在轴向位置处的径向换热对象的几何和热物性信息。节点轴向序号自井口向井底递增，井口节点轴向序号为O。(I)轴向网格划分轴向网格划分时，将求解区域依据隔水管深度、井身结构和钻具结构进行轴向分段，分段原则为:每一段内涉及到的所有传热对象的几何尺寸只有一种，分段信息记录在Sections数组中。然后自上至下对每一段再根据段长进行轴向网格划分，每一段内的网格控制体积相等。算法中使用NodeZ和NodeDZ两个数组记录网格中心轴向坐标和网格轴向长度。为节省程序内存，使用NodeInSection数组记录网格所在的段索引，大大减小了节点径向信息存储需要开支的存储空间。(2)径向信息记录从上到下通过扫描钻具结构和井身结构数据，分析每一段涉及到的径向传热对象几何信息和介质种类，记录在Annulus 二维数组中。假如某段有两层套管，套管外均为水泥环，则井筒内径向传热对象有:最内层为钻柱内钻井液，往外依次为钻柱、钻柱外环空钻井液、最内层套管1、套管I外水泥环、外层套管I1、套管II外水泥环,需要记录的信息有:钻柱内钻井液的介质类型(钻井液)、钻柱内钻井液的径向尺寸、钻柱管体介质类型(钢材)、钻柱管体径向尺寸、环空钻井液介质类型(钻井液)、环空钻井液径向尺寸、套管I介质(钢材)、套管I径向尺寸、套管I环空介质类型(水泥石)、套管I环空径向几何尺寸、套管II介质类型(钢材)、套管II径向几何尺寸、套管II环空介质类型(水泥石)、套管II环空径向几何尺寸。2)确定网格垂向坐标根据井眼轨迹确定网格中心的垂直坐标和网格垂直长度。3)应用初始条件海水区(泥线以上)钻柱内和环空所有节点的初始温度为节点垂直坐标对应深度处海水温度；地层区(泥线以下)钻柱内和环空所有节点的初始温度为网格中心垂直坐标对应垂直深度处地层原始温度；节点初始压力为气液两相钻井液静止状态下对应深度处的静液柱压力。4)自上而下计算钻柱内气液两相流钻井液节点温度和压力数据。每个节点的温度和压力都需要迭代计算至获得收敛解，具体计算步骤如下:①用变量TOld和POld记录上次迭代结束时该节点的温度和压力数据；②假定该节点的温度和压力等于上部相邻节点的温度和压力；③取该节点与上部相邻节点压力的平均值为网格单元平均压力；④确定环空气液两相流钻井液与钻柱外壁的强迫对流换热系数W2 ；⑤确定气相和液相在网格中心温度和平均压力下的热物性参数:密度、定压比热、热导率、动力粘度等；气相要根据具体气体类型是空气还是氮气，选择合适的高温高压热物性参数，液相则需根据钻井液类型选择合适的热物性计算方法。⑥按式(I)确定钻柱内气液两相流钻井液的压降梯度(dp/dZh ；
权利要求
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和压力差值H)iffMaX，以判断全局迭代计算是否满足结束条件。
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函数计算气液两相流摩阻压降， HeatResistancel:计算钻柱内钻井液到环空钻井液的换热热阻,实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数；HeatResistance2:计算环空钻井液到地层的换热热阻，实现函数功能需要调用HTPPipe和HTPAnnulus函数计算气液两相流在管流和环空流两种工况下的强迫对流换热系数； TInDrillStem函数根据钻柱内气液两相钻井液能量守恒方程迭代计算钻柱内气液两相钻井液节点温度； Tannulus:根据环空气液两相钻井液能量守恒方程迭代计算环空气液两相钻井液节点温度; ThermalPhysics:计算气相在给定温度和压力下的密度和比j:含,该函数功能需要根据气体类型调用空气； EOS:根据气体类型调用相应的状态方程，根据温度和比体积计算气相的压力和比焓； RKS:根据RKS模型计算给定温度和比体积时氮气的压力和比焓； AirTP:计算空气在给定温度和比体积时氮气的压力和比焓，计算比焓时需要调用AirAOT> AirArT> AirArDen函数计算空气状态方程的相关偏导数； AirAOT:计算空气的理想状态Helmholtz能对对比温度倒数的导数； AirArT:计算空气余能对对比温度倒数的导数； AirArDen:计算空气余能对对比密度的导数； AirConductivity:计算空气给定温度和压力下的热导率，实现函数功能需要调用ThermalPhysics函数计算空 气的真实密度； AirViscosity:计算空气给定温度和压力下的动力粘度，实现函数功能需要调用ThermalPhysics函数计算空气的真实密度； CpAir:计算空气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算空气的高温高压比j:含； CpNitrogen:计算氮气给定温度和压力下的定压比热，实现函数功能需要调用ThermalPhysics函数计算氮气的高温高压比j:含， NitrogenConductivity:计算氮气给定温度和压力下的热导率,实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度； NitrogenViscosity:计算氮气给定温度和压力下的动力粘度,实现函数功能需要调用ThermalPhysics函数计算氮气的真实密度。
全文摘要
本发明公开了一种深水气液两相流循环温度压力耦合计算方法，包括以下步骤1)计算网格垂向坐标；2)应用初始条件；3)从初始值出发，按照先钻柱内钻井液后环空钻井液的顺序迭代计算钻柱内和环空钻井液节点温度和压力数据，直至温度和压力都达到收敛条件，迭代结束，最后一次迭代计算结果为最终深海气液两相流井筒温度和压力模拟结果，保存并输出，该迭代命名为全局迭代。本发明的计算方法提高了计算精度。
文档编号G06F17/50GK103226641SQ20131016993
公开日2013年7月31日 申请日期2013年5月10日 优先权日2013年5月10日
发明者宋洵成, 管志川 申请人:中国石油大学(华东)