本发明是关于一种基于井筒模型的天然气井口流量计算方法，属于流体计算领域。
背景技术：
随着油气田开发不断向深水及沙漠等无人区迈进，对于油气集输系统而言新的问题不断出现，在传统的石油天然气工业中，为了简单高效实时地测量油井的产油、产气量和含水量，通常在每口井上安装一台多相流量计。
然而，在水下及无人值守的生产工艺诞生之后，传统的计量技术面临诸多新的问题，如果选择在水下或沙漠深处的生产系统中安装传统的多相流量计进行计量，该多相流量计不仅存在费用昂贵的问题，在深海或沙漠深处的条件下还存在日常的标定及维护困难的问题。
技术实现要素：
针对上述问题，本发明的目的是提供一种费用低廉且能够方便在深海或沙漠深处的日常标定及维护的基于井筒模型的天然气井口流量计算方法。
为实现上述目的，本发明采取以下技术方案：一种基于井筒模型的天然气井口流量计算方法，其特征在于，包括以下步骤：步骤1)：根据实际天然气生产井井筒的结构，建立天然气组份模型，并计算得到各压力和温度条件下生产井井口的产出物在流动过程中的物性参数；步骤2)：将井筒按照倾角和长度划分网格形成若干管段，并假设每一管段各截面处的质量流量相等；步骤3)：设定井筒出口压力和井筒出口温度的误差范围；步骤4)：根据实际测量的井筒入口温度、井筒入口压力、井筒出口温度和井筒出口压力，预测井筒的质量流量范围，并将该质量流量范围中的最低质量流量和最高质量流量的均值设定为初始的质量流量迭代值；步骤5)：设定每一管段的下游压力，根据物性参数、实际测量的井筒入口压力和井筒入口温度以及预设的管段下游压力、质量流量的迭代值和井筒出口压力的误差范围，自井筒上游至下游逐段推进，对已建立的井筒的热力模型和水力模型进行求解，迭代计算得到井筒出口温度和井筒平均持液率；步骤6)：根据预设的井筒出口温度误差范围，对计算的井筒出口温度与实际测量的井筒出口温度进行误差检验，若满足误差检验，则确定初始的质量流量迭代值为所求的气液混合质量流量，进入步骤7)；若不满足误差检验，则重新设定质量流量的迭代值，进入所述步骤5)；步骤7)：根据天然气组份模型以及计算的井筒平均持液率，将得到的气液混合质量流量进行拆分得到天然气生产井井口产出物中油、气和水各自的质量流量；步骤8)：将得到的油、气和水的质量流量分别转化为标况下的体积流量，完成天然气生产井井口产出物的流量计算。
进一步，在进行所述步骤1)～8)的天然气井口流量计算前，建立井筒的热力模型和水力模型，其中，水力模型包括持液率模型和压降模型，具体过程为：
A)建立井筒的热力模型
对于每一管段，将气相和液相视为整体，通过混合能量方程描述热力过程，并在混合能量方程的基础上进一步简化得到温降的简化解析解，在稳态流动的假设下不考虑摩擦力做功，建立稳态的混合能量方程：
式中，hm为气液混合物焓值；T为温度；P为压力；x为管段轴向；q为漏热量；g为重力加速度；θ为管段与水平面的夹角，在混合能量方程中：
式中，U为总传热系数；D0为管段外径；Te为环境温度；Wm为井筒的气液混合质量流量即质量流量的迭代值；气液混合物的定压质量比热容cpm和气液混合物的Joule-Thomson效应系数αhm分别为：
热力学参量的循环关系为：
将公式(2)～(5)代入混合能量方程(1)中得到：
将公式(6)整理后得到热力模型：
式中，参数
B)建立井筒的水力模型
a)建立井筒的持液率模型
假设气液相界面为平面，且气液相处于受力平衡状态，建立动量方程组：
式中，Ag为管段中气相所占的流通截面积；Al为管段中液相所占的流通截面积；τwg为管段中气相与壁面的剪切力；τwl为管段中液相与壁面的剪切力；τi为管段中气液相界面的剪切力；Sg为管段中气相的湿周长度；Si管段中气液相界面的湿周长度；Sl为管段中液相的湿周长度；ρg为管段中气相的密度；ρl为管段中液相的密度；
将动量方程组消去压力梯度整理后得到联合动量方程即持液率模型：
b)建立井筒的压降模型
将动量方程组中的两动量方程相加得到井筒的压降模型：
进一步，所述步骤4)中井筒的质量流量范围采用IPR方程，根据实际测量的井筒入口温度、井筒入口压力、井筒出口温度和井筒出口压力得到。
进一步，所述步骤5)中设定每一管段的下游压力，根据物性参数、实际测量的井筒入口压力和井筒入口温度以及预设的管段下游压力、质量流量的迭代值和井筒出口压力的误差范围，自井筒上游至下游逐段推进，对已建立的井筒的热力模型和水力模型进行求解，迭代计算得到井筒出口温度和井筒平均持液率，具体过程为：步骤5.1)：设定管段的下游压力，根据物性参数、实际测量的井筒入口压力和井筒入口温度以及预设的管段下游压力和质量流量的迭代值，对已建立的井筒的热力模型进行求解，计算得到管段的下游温度Ti：
式中，Ti-1为该管段的上游温度即上一管段的下游温度，且井筒入口处管段的上游温度和上游压力即为实际测量的井筒入口温度和井筒入口温度；步骤5.2)：根据管段的分层流几何关系，对已建立的井筒的持液率模型进行求解，计算得到该管段的当量液面高度和平均持液率；步骤5.3)：根据管段的分层流几何关系、摩阻封闭关系和得到的当量液面高度，对已建立的井筒的压降模型进行求解，计算得到该管段的下游压力；步骤5.4)：重复所述步骤5.1)～5.3)，自井筒上游至下游逐段推进，迭代求解得到井筒中每一管段的下游温度、下游压力和平均持液率，其中，最后一次预设的管段下游压力为实际测量的井筒出口压力；步骤5.5)：根据预设的井筒出口压力误差范围，对最后一次迭代得到的下游压力与最后一次预设的管段下游压力进行误差检验，若满足误差检验，则将最后一次迭代得到的下游温度、下游压力以及平均持液率视为井筒出口温度、井筒出口压力和井筒平均持液率，进入所述步骤6)；若不满足误差检验，则重新设定管段的下游压力，进入所述步骤5.1)。
进一步，所述步骤5.2)中管段的分层流几何关系包括：管段中液相的体积含率与当量液面高度的几何关系：
式中，γ为管段中分层流液相所对应的圆心角；hl为管段中液相的液位高度；d为管段内径；αl为管段中液相的体积含率即管段的平均持液率；管段中气相和液相的湿周长度与当量液面高度的几何关系：
管段中气相和液相所占流通截面积与当量液面高度的几何关系：
式中，A为管段中气相和液相的总流通截面积；管段中气相和液相水力直径的几何关系：
式中，dl为管段中液相的水力直径；dg为管段中气相的水力直径。
进一步，所述步骤5.3)中管段的摩阻封闭关系为：通过范宁公式计算气液相、管壁以及气液相之间的范宁摩阻系数：
式中，fg为管段中气相的范宁摩阻系数；fl为管段中液相的范宁摩阻系数；fi为管段中气液相界面的范宁摩阻系数；vg为管段中气相的速度；vl为管段中液相的速度；管段中气液相界面的范宁摩阻系数fi在层流时采用Hagen-Poiseill相关式进行计算，即：
式中，Rei为管段中气液相界面的雷诺数，且Rei＞2300；管段中气液相界面的范宁摩阻系数fi在湍流时采用Colebrook-White相关式进行计算，即：
式中，di为管段中气液相界面的水力直径，且Rei≤2300；对于管段中气液相界面的范宁摩阻系数fi，认为其与管段中气相的范宁摩阻系数相同fg。
进一步，所述步骤6)中根据预设的井筒出口温度误差范围，对计算的井筒出口温度与实际测量的井筒出口温度进行误差检验，若满足误差检验，则确定初始的质量流量迭代值为所求的气液混合质量流量，进入步骤7)；若不满足误差检验，则进入所述步骤5)重新设定质量流量的迭代值，具体过程为：若计算的井筒出口温度和实际测量的井筒出口温度的相对误差在±1×10-4之间，则初始的质量流量迭代值满足误差检验的精度要求，为所求的气液混合质量流量，进入所述步骤7)；若计算的井筒出口温度和实际测量的井筒出口温度的相对误差大于1×10-4，则预测的质量流量偏大，采用二分法将最低质量流量和初始的质量流量迭代值的平均值作为质量流量的迭代值，进入所述步骤5)；若计算的井筒出口温度和实际测量的井筒出口温度的相对误差小于-1×10-4，则预测的质量流量偏小，采用二分法将最高质量流量和初始的质量流量迭代值的平均值作为质量流量的迭代值，进入所述步骤5)。
本发明由于采取以上技术方案，其具有以下优点：本发明的天然气井口流量计算方法只需要实际测量井筒入口温度、井筒入口压力、井筒出口温度和井筒出口压力，其他所需参数均可通过一系列计算得出，不需要额外的实体流量计量仪表就能够计算出天然气生产井井口产出物各分量的流量，可对单个或多个天然气井实现油、气。水三相的分相计量，克服了传统的分离计量和多相流量计计量方法既难于应用于水下和沙漠深处的生产系统，也难于维护和标定的难题，可以广泛应用于海上或沙漠天然气开采行业领域中。
附图说明
图1是本发明中天然气组份模型的结构示意图，其中，图1(a)为井筒内分层流动的物理模型主视图，图1(b)为图1(a)的截面图；
图2是本发明天然气井口流量计算方法的流程图；
图3是本发明中管段的分层流几何关系的原理示意图。
具体实施方式
以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。
在进行计算之前，需要建立井筒的热力模型和水力模型，其中，水力模型包括持液率模型和压降模型,具体为：
1)建立井筒的热力模型
对于每一管段，将气相和液相视为整体，通过混合能量方程描述热力过程，并在混合能量方程的基础上进一步简化得到温降的简化解析解。在稳态流动的假设下不考虑摩擦力做功，建立稳态的混合能量方程：
式中，hm为气液混合物焓值；T为温度；P为压力；x为管段轴向；q为漏热量；g为重力加速度；θ为管段与水平面的夹角。
在混合能量方程(1)中：
式中，U为总传热系数；D0为管段外径；Te为环境温度；Wm为井筒的气液混合质量流量。气液混合物的定压质量比热容cpm和气液混合物的Joule-Thomson(焦耳-汤姆孙)效应系数αhm分别为：
热力学参量的循环关系为：
将公式(2)～(5)代入混合能量方程(1)中得到：
将上述公式(6)整理后得到热力模型：
式中，参数
2)建立井筒的水力模型
2.1)建立井筒的持液率模型
如图1所示，处于分层流时每一管段的气相和液相均有各自的流通通道，气相和液相间的相互作用并不强烈，假设气液相界面为平面，且气液相处于受力平衡状态，建立动量方程组：
式中，Ag为管段中气相所占的流通截面积；Al为管段中液相所占的流通截面积；τwg为管段中气相与壁面的剪切力；τwl为管段中液相与壁面的剪切力；τi为管段中气液相界面的剪切力；Sg为管段中气相的湿周长度；Si管段中气液相界面的湿周长度；Sl为管段中液相的湿周长度；ρg为管段中气相的密度；ρl为管段中液相的密度。
将动量方程组(8)消去压力梯度整理后得到联合动量方程即持液率模型：
2.2)建立井筒的压降模型
将动量方程组(8)中的两动量方程相加得到井筒的压降模型：
如图2所示，基于建立的井筒的热力模型和水力模型，本发明提供的基于井筒模型的天然气井口流量计算方法，包括以下步骤：
1)根据实际天然气生产井井筒的结构,建立天然气组份模型，并计算得到各压力和温度条件下生产井井口的产出物(油、气和水)在流动过程中的物性参数，其中，物性参数包括比热容、传热系数和密度等，建立和求解天然气组份模型可以采用现有技术公开的天然气组份建模及求解方法，具体过程在此不做赘述。
2)将井筒按照倾角和长度划分网格形成若干管段，并假设井筒中流体的流动过程是稳态流动过程，即每一管段各截面处的质量流量相等，其中，网格的精度可以根据实际情况进行设定。
3)设定井筒出口压力和井筒出口温度的误差范围，均可以根据实际情况所需进行设定。
4)根据与气田生产数据采集系统(DCS)建立通信连接获取的常规生产中实际测量的井筒入口温度、井筒入口压力、井筒出口温度和井筒出口压力，采用IPR(Inflow Performance Relationship油/或气井流入动态特性)方程预测井筒的质量流量范围，并将该质量流量范围中的最低质量流量和最高质量流量的均值设定为初始的质量流量迭代值，质量流量的迭代值即为井筒的气液混合质量流量Wm，其中，采用IPR方程预测井筒的质量流量范围为现有技术，具体过程在此不做赘述。
5)设定每一管段的下游压力，根据物性参数、实际测量的井筒入口压力和井筒入口温度以及预设的管段下游压力、质量流量的迭代值和井筒出口压力的误差范围，自井筒上游至下游逐段推进，对已建立的井筒的热力模型和水力模型进行求解，迭代计算得到井筒出口温度和井筒平均持液率，具体为：
5.1)设定管段的下游压力，根据物性参数、实际测量的井筒入口压力和井筒入口温度以及预设的管段下游压力和质量流量的迭代值，对已建立的井筒的热力模型进行求解，计算得到管段的下游温度，其中，管段的下游压力可以根据实际情况进行设定。
对于端部节点为(i-1，i)(端部节点中的i表示井筒管段的编号)、长为Δx、与水平面夹角为θ的管段，认为其参数U、cpm、αhm、Te和均为常量，对热力模型(7)在Δx进行积分，求解得到该管段的下游温度Ti：
式中，Ti-1为该管段的上游温度即上一管段的下游温度，且井筒入口处管段的上游温度和上游压力即为实际测量的井筒入口温度和井筒入口温度。
5.2)根据管段的分层流几何关系，对已建立的井筒的持液率模型进行求解，计算得到该管段的当量液面高度和平均持液率。
如图3所示，根据分层流几何关系，持液率模型(9)中的所有变量均为当量液面高度的函数，因此持液率模型(9)实际上为管段的液面高度隐式方程，即持液率方程。
通过二分法求解该持液率方程即能够得到无量纲的当量液面高度并根据下述管段的分层流几何关系，求解得到管段的平均持液率。
管段的分层流几何关系包括：
管段中液相的体积含率与当量液面高度的几何关系：
式中，γ为管段中分层流液相所对应的圆心角；hl为管段中液相的液位高度；d为管段内径；αl为管段中液相的体积含率即管段的平均持液率。
管段中气相和液相的湿周长度与当量液面高度的几何关系：
管段中气相和液相所占流通截面积与当量液面高度的几何关系：
式中，A为管段中气相和液相的总流通截面积。
管段中气相和液相水力直径的几何关系：
式中，dl为管段中液相的水力直径；dg为管段中气相的水力直径。
5.3)根据管段的分层流几何关系、摩阻封闭关系和得到的当量液面高度，对已建立的井筒的压降模型进行求解，计算得到该管段的下游压力。
气液相、管壁以及气液相之间的摩阻封闭关系(通过范宁公式得到)：
式中，fg为管段中气相的范宁摩阻系数；fl为管段中液相的范宁摩阻系数；fi为管段中气液相界面的范宁摩阻系数；vg为管段中气相的速度；vl为管段中液相的速度。
其中，管段中气液相界面的范宁摩阻系数fi在层流时采用Hagen-Poiseill(哈根-泊肃叶)相关式进行计算，即：
式中，Rei为管段中气液相界面的雷诺数，且Rei＞2300。
管段中气液相界面的范宁摩阻系数fi在湍流时采用Colebrook-White(柯尔勃洛克-怀特)相关式进行计算，即：
式中，di为管段中气液相界面的水力直径，且Rei≤2300。
对于管段中气液相界面的范宁摩阻系数fi，认为其与管段中气相的范宁摩阻系数相同fg。
求解管段的下游压力时，通过持液率模型(9)得到管段的当量液面高度后，再根据管段的分层流几何关系和摩阻封闭关系计算得到压降模型(10)右侧的变量值，进而得到该管段的下游压力。
5.4)重复步骤5.1)～5.3)，自井筒上游至下游逐段推进，迭代求解得到井筒中每一管段的下游温度、下游压力和平均持液率，其中，最后一次预设的管段下游压力为实际测量的井筒出口压力。
5.5)根据预设的井筒出口压力误差范围，对最后一次迭代得到的下游压力与最后一次预设的管段下游压力进行误差检验，若满足误差检验，则将最后一次迭代得到的下游温度、下游压力以及平均持液率视为井筒出口温度、井筒出口压力和井筒平均持液率，进入步骤6)；若不满足误差检验，则重新设定管段的下游压力，进入步骤5.1)。
例如：若最后一次迭代得到的下游压力与最后一次预设的管段下游压力的绝对误差在±5000kPa之间，则最后一次迭代得到的下游压力满足误差检验的精度要求，最后一次迭代得到的下游温度、下游压力以及平均持液率为所求的井筒出口温度、井筒出口压力和井筒平均持液率，进入步骤6)。若最后一次迭代得到的下游压力与最后一次预设的管段下游压力的绝对误差大于5000kPa或小于-5000kPa，则最后一次迭代得到的下游压力不满足误差检验的精度要求，需要重新设定管段的下游压力，进入步骤5.1)。
6)根据预设的井筒出口温度误差范围，对计算的井筒出口温度与实际测量的井筒出口温度进行误差检验，若满足误差检验，则确定初始的质量流量迭代值为所求的气液混合质量流量，进入步骤7)；若不满足误差检验，则采用二分法重新设定质量流量的迭代值，进入步骤5)，具体为：
若计算的井筒出口温度和实际测量的井筒出口温度的相对误差在±1×10-4之间，则初始的质量流量迭代值满足误差检验的精度要求，为所求的气液混合质量流量，进入步骤7)。若计算的井筒出口温度和实际测量的井筒出口温度的相对误差大于1×10-4，则说明预测的质量流量偏大，采用二分法将最低质量流量和初始的质量流量迭代值的平均值作为质量流量的迭代值，进入步骤5)；若计算的井筒出口温度和实际测量的井筒出口温度的相对误差小于-1×10-4，则说明预测的质量流量偏小，采用二分法将最高质量流量和初始的质量流量迭代值的平均值作为质量流量的迭代值，进入步骤5)。
7)根据天然气组份模型以及井筒平均持液率，将得到的气液混合质量流量在实际测量的井筒入口压力和温度条件下，通过闪蒸计算拆分得到天然气生产井井口产出物中油、气和水各自的质量流量，其中，闪蒸计算可以采用现有技术公开的方法，具体过程在此不做赘述。
8)将得到的油、气和水的质量流量分别转化为标况下的体积流量，完成天然气生产井井口产出物的流量计算。
本发明在实际使用时，可以将本发明的天然气井口流量计算方法视为一个抽象函数F，当井筒的物理参数(井筒路由、组分、传热等)已知，井筒入口温度和井筒入口压力作为固定的输入参数，则该抽象函数F成为井筒气液混合质量流量的单值函数，即输入一个井筒的气液混合质量流量便能得到所对应的井筒入口、出口的温度差和压力差，因此，存在以下方程：
F(P入口实测，T入口实测，Wm)＝ΔT＝T入口实测-T出口实测 (19)
式中，P入口实测为实际测量的井筒入口压力，T入口实测为实际测量的井筒入口温度，T出口实测为实际测量的井筒出口温度。采用二分法求解上述非线性方程(19)，即可以得到井筒的气液混合质量流量。关于二分法迭代的上下限，可以通过IPR方程得到该工况下的近似产能(即井筒的气液混合质量流量)，并将该近似产能值的4倍作为迭代的上限，该近似产能值的四分之一作为迭代的下限。以本发明的天然气井口流量计算方法为核心编制的计算机程序可以实现在线计算天然气井的总流量、气相、液相的分相流量等生产数据，可以在一些情况下代替传统方法所需的实体流量计量仪表。
上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。