目前，国内各油田采用的油井产量计量方法主要有：玻璃管量油孔板测 气、翻斗量油孔板测气、两相分离密度法和三相分离计量方法等。
1.原油的测量
采用玻璃管液面计量油：
玻璃管量油是国内各油田普遍采用的传统方法，约占油井总数的90％以 上。该方法的有益之处在于装备简单、投资少；但是，由于采用间歇量油的 方式来折算产量，导致原油系统误差较大，约为10％～20％。
采用电报量油：
在玻璃管液面计量油的基础上，在规定的量油高度H上、下各安装一个电 极，当水上升到下电极时，计时电表接通开始计时；水上升到上电极时，电 表切断停止走动，记录水上升H高的时间t，则可按照玻璃管液面计量油的方 法计算出油井的产量。该方法的主要缺点是操作难度较大，不能连续量油， 因为需要人工导流程，仍然需要每次量完后用其它井的天然气将计量罐中的 液面压下来才能进行下一次计量。
采用翻斗量油：
翻斗量油装置主要由量油器、计数器等组成。一个斗装满时翻到排油， 另一个斗装油，这样反复循环来累积油量。这种量油装置结构简单，具有一 定计量精度。该方法的主要缺点在于，设备投资较大，不能应用普通的计量 罐，而且经常出现漏量或倒不干净的情况，对于产量很低的稠油适应性稍好 些。
2.综合测量的计量装置
随着技术的发展以及油田管理和降低工人劳动强度、提高生产效率的需 要，相继出现了许多可以对油井油、气、水产量进行综合测量的计量装置。
采用三相分离计量：
三相分离计量是把油、气、水分离后分别计量，分离后原油含水较低(一 般在30％以下)，原油测量误差降低，不受油井含水率的影响。但是，要想把 特高含水原油分离成低含水原油并进行计量，工艺技术十分复杂，而且数量 很大的游离水经常携带一部分原油，造成很大误差，且所需的设备仪表多、 投资大、管理操作难度大、维修费用高。
采用两相分离计量：
两相分离计量是将油井采出液分离成液体和气体，然后分别对其进行计 量。两相分离计量设备主要由两相分离器、气体流量计、液体流量计、含水 分析仪等组成。其中，气体流量计和液体流量计计量油井的产气量和产液量， 含水分析仪测量分离出液体的含水率，由此计算出油井的油、气、水产量。
另一种形式的计量设备由两相分离器、质量流量计和气体流量计组成。 质量流量计测量分离出的液量，并计算出其中的含水率，从而测量出油井的 油、气、水产量。这种计算装置投资较少、操作简便，在我国油田中获得了 较多的应用。该方法的主要缺点是需要较复杂的地面管线和流程，地面管线 和建站投资大。
采用不分离计量：
不分离计量是不分离油井采出液，将文丘里管、密度计或不同的流量传 感器结合起来计量气体和液体的流量，液体部分用双γ射线密度计、电容、 微波水含量监控仪来确定油和水的含量，从而计算出油、气、水各自的产量。
油气水不分离计量技术在占地面积等方面有很大的优越性。但是，油井 采出液中的油、气、水等组分一般不是均匀混合的，并以不同的速度流动， 还可能相互作用形成蜡和氢氧化物，并且引起难以预料的复杂流态。因此， 开发具有广泛适用范围的流量计具有很大的难度，并且成功投入使用的还很少。
并且，在上述方法中，均需设置油井计量站，并且通过人工定期操作录 取数据，这样不可避免地存在非连续性和需建立相关的地面计量流程及装置 等问题，同时增加了成本。

鉴于现有技术中存在的问题，本发明的目的在于提供一种油井产液量计 量、工况分析优化方法及其系统。通过本发明，可实现在油井无计量站的情 况下，实现远程油井工况数据的自动录取、在无人值守情况下及时掌握油井 的动态变化。
本发明提供一种油井产液量计量、工况分析优化方法，该方法包括步骤：
获取油井抽油机上的传感器传送的工况数据，并通过无线通信网络将该 工况数据传送至工况采集及监视单元；
工况采集及监视单元接收所述工况数据，并将该工况数据进行处理后传 送至产液量计量单元，并监视油井的运行状态；
产液量计量单元接收到所述工况数据后，根据该工况数据和数据库中存 储的油井基础数据计算产液量。
根据该方法，还包括步骤：
对工况数据和计算出的产液量进行数据分析；
根据数据分析结果进行工况分析；
根据数据分析和工况分析结果进行优化设计。
本发明还提供一种油井产液量计量、工况分析优化系统，该系统至少包 括：数据采集控制器、工况采集及监视单元、产液量计量单元和存储单元； 其中，数据采集控制器，安装于油井中并与设置在油井抽油机上的传感器相 连接，用于采集所述传感器的工况数据并通过无线通信网络将该工况数据传 送至工况采集及监视单元，并对油井进行控制；
工况采集及监视单元，通过无线通信网络与数据采集控制器进行信息交 互，接收所述数据采集控制器发送的工况数据，并将手动或自动设置的指令 发送至数据采集控制器，并监视油井的运行状态；
产液量计量单元，与所述工况采集及监视单元和存储器连接，接收所述 工况采集及监视单元发送的工况数据，并根据该工况数据和储存在存储单元 中的油井基础数据，并且应用产液量计算数学模型计算产液量；并将采集到 的工况数据送入存储器进行存储；
存储单元，与所述产液量计量单元连接，用于储存油井基础数据供产液 量计量单元使用；接收产液量计量单元传送的产液量计量结果并储存，并接 收采集的工况数据并储存。
还包括数据分析单元，与所述产液量计量单元和存储单元连接，用于对 相关数据进行分析并将数据分析结果发送至存储器进行数据储存；其中，所 述相关数据至少包括：产液量、压力、载荷、泵效、系统效率、区块中单井 运行合理性统计。
还包括工况分析单元，与所述数据分析单元连接，根据数据分析单元的 分析结果对油井工况进行分析。
还包括优化设计单元，与所述数据分析单元和工况分析单元连接，根据 数据分析结果和工况分析结果对油井进行优化设计。
本发明的有益效果在于，可在油井无计量站情况下实现远程油井工况数 据的自动录取、在无人值守情况下能及时掌握油井的动态变化；
低成本投入、高可靠性和易维护、可拓展；
以油井单井产量计量为核心，并辅以对油井进行工况监控，对采集的数 据进行分析诊断优化，替代或简化计量流程，以降低产能建设投入和运行成 本的同时，还实现提高油井系统效率的目的。
附图说明
图1为本发明实施例1的油井产液量计量、工况分析优化系统构成示意图；
图2为游梁式抽油机传感器安装及数据采集控制器安装示意图；
图3为螺杆泵抽油机传感器安装及数据采集控制器安装示意图；
图4为电泵抽油机传感器安装及数据采集控制器安装示意图；
图5为游梁机有杆泵井产液量计算流程图；
图6为自喷井产液量计算流程图；
图7为螺杆泵井产液量计算流程图；
图8A和图8B为梁的动力学方程理论示意图；
图9为电潜泵井产液量计算流程图；
图10为求粘度准数系数CNt的相关关系图；
图11为滞留系数相关关系图；
图12为再次校正系数相关关系图。

以下结合附图对本发明进行详细说明。
实施例1
本发明提供一种油井产液量计量、工况分析优化系统。如图1所示，该系 统至少包括：数据采集控制器(RTU)101、工况采集及监视单元102和产液量 计量单元103；其中，
数据采集控制器RTU101，安装于油井中并与设置在油井抽油机上的传感 器相连接，用于采集所述传感器的工况数据并通过无线通信网络将该工况数 据传送至工况采集及监视单元102，并对油井进行控制；其中，对油井进行控 制是指对油井进行开井、停井、发出开停井声音警报等智能控制；
工况采集及监视单元102，通过无线通信网络与数据采集控制器101进行 信息交互，接收所述数据采集控制器101发送的工况数据，并将手动或自动设 置的指令发送至数据采集控制器，并监视油井的运行状态；其中，工况采集 及监视单元102与数据采集控制器101的信息交互具体是指：工况采集及监视 单元102不但接收数据，而且可以与数据采集控制器RTU101进行信息交互，将 中心控制室的人工输入的或自动设置的开井、停井、开停井警报等指令传递 给数据采集控制器RTU101，由数据采集控制器RTU101控制执行开井、停井、 开停井警报等工作；
产液量计量单元103，与所述工况采集及监视单元102和存储单元107连 接，接收所述工况采集及监视单元102发送的工况数据，并根据该工况数据和 储存在存储单元107中的油井基础数据，应用油井产液量计算数学模型计算产 液量；并且将采集的工况数据存入存储单元107中；
存储单元107，与所述产液量计量单元103连接，用于储存油井基础数据 供产液量计量单元103使用；接收产液量计量单元103传送的产液量计量结果 并存储；并且接收数据分析单元104传送的数据分析结果。
其中，无线通信方式可采用GSM、GPRS或CDMA通信方式。
本实施例中，数据采集控制器101采用自主开发的YDSW远程数据采集控制 器RTU，其为高精度的数据采集器，可对油井进行集中控制和自动管理，可无 线监测控制自喷井、游梁机有杆泵井、电潜泵井、螺杆泵井、注水井生产状 况。
如图2所示，数据采集控制器101获取游梁抽油机上的载荷(负荷)传感 器201、位移传感器202、压力(油压)传感器203、温度传感器204、曲柄脱 扣传感器205等传送的数据；此外还包括通信单元206、监控箱207、配电柜208 (内有高压转换单元、电压/电流/CAN(控制器局域网)单元)、井场照明灯 209和红外探测器210。
如图3所示，数据采集控制器101获取螺杆泵抽油机上的扭矩载荷一体化 传感器和转速传感器301、温度传感器302、压力(油压)传感器等传送的数 据。此外还包括电控箱304，箱内有通信单元、电压/电流/CAN单元、高压转 换单元。
如图4所示，数据采集控制器101获取电潜泵抽油机上的油嘴前压力传感 器401，油嘴后压力传感器402、和套压变送器403(图中未示出)传送的数据； 此外还包括电控箱404，箱内有通信单元、电压/电流/CAN单元、高压转换单 元。
对于自喷井，数据采集控制器101可获取安装在自喷井上的油压、回压变 送器和套压变送器发送的数据。
所述工况采集及监视单元102可为一台或多台计算机，采集所安装传感器 的油井的各种工况数据，并监视油井运行状态。该工况采集及监视单元102通 过无线通信方式接收所述数据采集控制器101传送的数据；也将中心控制室的 人工输入的或自动设置的开井、停井、开停井警报等指令传递给数据采集控 制器RTU101，由数据采集控制器RTU101控制执行开井、停井、开停井警报等 工作。
另外，该工况采集及监视单元102还可通过一台服务器与数据采集控制器 RTU101交互信息，此服务器可为油田公司信息中心Internet主机。
本实施例中还包括数据分析单元104，与所述产液量计量单元103连接， 用于对相关数据进行分析。其中，相关数据是指：油井的产液量、压力、载 荷、泵效、系统效率、区块单井运行合理性统计等指标。
本实施例中，所述存储单元107包括：第一存储单元，用于存储油井基础 数据；第二存储单元，用于存储油井生产报表，该生产报表中包括产液量、 采集的压力、温度、转速、电参数等参数，还可包括应用数据分析单元104计 算的数据分析结果，如泵效组成、系统效率、局部损耗、耗电量、耗电成本、 油井诊断结果等技术指标。该存储单元107可采用数据库服务器实现。本实施 例中可为油田公司信息中心数据库服务器。
本实施例中，产液量计量单元103、数据分析单元104可在服务器上实现， 产液量计量单元103根据工况数据和储存在存储单元107中的油井基础数据， 应用自喷井、电潜泵井、螺杆泵井、游梁机有杆泵井产液量计算数学模型， 计算各类油井的单井产液量，并在数据分析单元104中进行数据分析，然后并 将产液量、采集的工况数据和数据分析结果一同传入存储单元107中的生产报 表，即存储单元107中的第二存储单元中进行储存。抽油机井、自喷井、电潜 泵井、螺杆泵井产液量计算的流程图如图5、图6、图7、图9所示。
还包括工况分析单元105，与所述数据分析单元104连接，根据数据分析 单元104的分析结果对油井工况进行分析。其中，对油井工况进行分析是指对 油井工况进行诊断，可以自动诊断常见的19种故障：连抽带喷、固定凡尔卡 死(不能打开)、泵严重磨损(不能关闭)、抽油杆断脱、气锁、完全液击、 气体影响、供液不足、柱塞脱出工作筒、固定凡尔漏失、游动凡尔漏失、液 体或机械摩阻、泵筒弯曲、泵上碰、泵下碰、卡泵、泵工作基本正常等油井 工作状况诊断分析。
还包括优化设计单元106，与所述工况分析单元105连接，根据工况分析 结果对油井进行优化设计。
本实施例中，还包括一个或多个用户终端108，与所述产液量计量单元103 连接并进行信息交互，对油井产液量信息进行维护、对产液量计量结果进行 查询以及依据数据分析与工况分析结论进行优化设计。
如图1所示，该系统还包括远程视频监控单元109，与所述工况采集及监 视单元102连接，通过无线通信网络对油井工况进行实时监控。即在油井、接 转站或联合站外安装至少一个云台和摄像机，可为全景低照度摄像机，对站 内全貌及油井进行监视。
该系统还包括网络浏览单元110，与所述远程视频监控单元109、工况采 集及监视单元102、产液量计量单元103、数据分析单元104、工况分析单元105、 优化设计单元106以及存储单元107连接，用于生产状况和相关数据进行实时 浏览、查询。本实施例中可通过IE浏览器和装有的视频播放软件，在油田信 息网上可随时浏览各油井的各个监控画面及实时生产数据，及产液量计算结 果，查询有关生产报表及分析结果。
上述系统中，油水井无线工况数据监控：数据采集控制器RTU，即无线油 水井工况智能远程监控器是对油井进行集中控制和自动管理的专用数据采集 器，其数据传输方式为先进可靠的无线数据传输。可无线监测控制自喷井、 游梁机井、螺杆泵井、电潜泵井、注水井生产工况。数据传输方式，可选用 GSM/GPRS/CDMA通信方式。
液量计算及分析优化：油水井液量计算、油水井优化设计、工况诊断、 工况宏观管理、系统效率分析与评价决策等于一体的综合软件，提供了一个 优秀的集成解决方案。
无线网络视频监视：引进国际领先的、基于MPEG4硬件压缩技术的网络 视频服务器，可以在64K-2M的低带宽占用情况下传输高清晰度的实时视频图 像，通常情况下200K-300K即可实现30帧/秒的全实时监视。产品采用TCP/IP 协议，可以基于局域网/广域网/Internet构建大规模的集中图像监控系统。
网络浏览：系统为浏览器/服务器(Browser/Server)构架。油田局域网 内的授权用户可进行实时浏览、查询工况数据。
实施例2
本发明还提供一种油井产液量计量、工况分析优化方法，该方法采用上 述计量系统完成，该方法包括步骤：
数据采集控制器101获取设置在油井抽油机上的传感器传送的工况数据， 并通过无线通信网络将该工况数据传送至工况采集及监视单元102；
工况采集及监视单元102接收所述工况数据并进行处理后，并将该工况数 据传送至产液量计量单元103，并监视油井的运行状态；其中，对所述工况数 据进行处理是指将工况数据进行加密处理并打包；油井的运行状态具体指： 是否开井、是否缺相、载荷位移是否采集得到等工作状况；另外，还可包括 通过数据采集控制器RTU101初算是否没有电流、电流是否缺相等指标；
产液量计量单元103接收到所述工况数据后，根据该工况数据和数据库中 存储的油井基础数据计算产液量；其中，油井基础数据为油密度、油粘度、 气油比、含水、动液面、泵径、泵深、杆柱组合等。
还包括步骤：
对工况数据和计算出的产液量进行数据分析，其中可包括产液量分析、 压力分析、载荷分析、泵效分析、区块中单井运行合理性统计分析等；根据 数据分析结果进行工况分析；根据工况分析结果进行优化设计。
还包括步骤：对采集的所述工况数据、计算出的产液量和数据分析结果 进行存储。
其中，所述根据工况数据和油井基础数据计算产液量，包括步骤：
根据工况数据和油井基础数据，应用产液量计算数学模型计算油井的单 井产液量；通过流量标定系数进行修正；
将修改正后的计算产液量作为油井计量产液量。
在产液量的计算中，由于产液量计算数学模型较为复杂，一些地质因素 无法考虑进去，有时不同的区块会或多或少存在系统误差，需要应用罐车计 量的标准产液量标定，得到流量标定系数，通过流量标定系数进行修正，剔 除系统误差，将修改正后的计算产液量作为油井计量产液量。
本实施例中，具体的标定系数计算可采用如下方法：
对于每口油井：
$K_i=\frac{q_g-q_y}{q_y}---(F0-1)$
整个区块：
$K=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{K}_i---(\mathrm{FD}-2)$
因此，通过上述标定系数K对油井产液量qg进行修正；
修正后的油井产液量：
Q＝Kqg                                            (F0-3)
式中：Q：油井产液量，K：修正系数，qg：功图法计算的产液量，qy： 实际产量。
上述标定系统解决办法适用于任何类型的油井。
下面分别以游梁机有杆泵井(抽油机井)、自喷井、螺杆泵井、电潜泵 井为例对本发明的产液量计量方法进行详细说明。
1.游梁机有杆泵井
如图5所示，为游梁机有杆泵井时计算产液量的流程图。其中，采用有杆 泵抽油时，根据地面示功图求得井下各级杆柱功图和泵功图；应用泵功图识 别技术计算产液量。
本模型中，所用到的工况数据为：瞬时产液量、累积产液量，泵功率、 光杆功率、泵效、漏失、气体影响泵效、有效功率，系统损耗、系统效率、 泵排出压力，平衡状况、平衡建议，电量、电费等，还有丰富的宏观控制图、 诊断结论等诸多指标。
所用到的油井基础数据为：示功图、三相电流I1，I2，I3、电压U、功率因素 cosφ、油压、套压等动态参数，油密度、油粘度、气油比、含水、动液面、 泵径、泵深、杆柱组合、生产气液比Rs等静态参数根据。
通过测试抽油机井所测的表示抽油杆柱最上端的载荷和位移之间关系的 地面示功图，应用示功图、油密度、油粘度、气油比、含水、泵径、泵深、 杆柱组合、生产气液比Rs等静态参数根据，应用杆柱、液柱和油管三维振动 数学模型求解，即有杆泵抽油系统模型，得到井下各级杆柱功图和泵功图； 应用泵功图识别技术计算油井产液量和各种工况指标；其中，先计算井下泵 的有效排量，然后再计算井口折算有效排量。
其理论模型示意函数如下式：
Q有杆泵＝kf(s，n，Dp，LP，GT，μ，Rs)
其中：Q有杆泵----有杆泵井的产液量，m3/d；s----冲程，m；n----冲次， r/min；DP----泵径，mm；LP----杆柱组合，即多极杆柱的每段的钢级、杆 径和杆长，m；GT----地面示功图数据，包括位移、载荷、冲次等数据；μ---- 井液粘度，mpa.s；Rs----生产气液比；k----流量标定系数，小数。
本实施例中，有杆泵抽油系统包括杆柱、液柱和油管三维振动，主要考 虑了在抽油过程中抽油杆、液柱和油管三个运动系统之间的相互关系，如图5 所示。它们受到多种力的作用，其中包括垂力(杆柱重、液柱重、油管重)、 惯性力、摩擦力(杆管之间摩擦力、杆液之间摩擦力、管液之间摩擦力、柱 塞副之间摩擦力等)、振动、沉没压力等。根据力的合成，可以得到描述油 管运动、抽油杆运动和液柱运动的偏微分方程组。
本部分计算模型较复杂，主要应用以下考虑了抽油杆、油管、液柱的三 维振动方程，其偏微分方程组如下：
杆柱的波动方程：

                           Pf(0，t)＝P0
按照如上初始条件和边界条件，利用数值方法求解这些偏微分方程组， 从而求得任意深度、任意时间的功图、压力分布、载荷、位移、速度等参数。 具体的偏微分方程组采用下述数值解法求解，但不限于下述解法。
本实施例中，可采用傅里叶级数解法，求得泵柱塞处的动载荷与柱塞位 移之间关系，具体方法如下。
用以截尾傅里叶级数表示的悬点动载荷函数D(t)及光杆位移函数U(t)作 为边界条件：
$D\left(t\right)=\frac{{\sigma }_0}{2}+\underset{n=1}{\overset{\stackrel{-}{n}}{\Sigma }}({\sigma }_n\cos \mathrm{\omega t}+{\tau }_n\sin \mathrm{n\omega t})---(F1-1)$
$U\left(t\right)=\frac{v_0}{2}+\underset{n=1}{\overset{\stackrel{-}{n}}{\Sigma }}(v_n\cos \mathrm{n\omega t}+{\delta }_n\sin \mathrm{n\omega t})---(F1-2)$
D(t)及U(t)的傅里叶系数σoσnτn及vovnδn可分别用下面的公式求得：
${\sigma }_n=\frac{\omega }{\pi }{\int }_0^{T}D\left(t\right)\cos \mathrm{n\omega tdt}(n=\mathrm{0,1,2},...\stackrel{-}{n})---(F1-3)$
${\tau }_n=\frac{\omega }{\pi }{\int }_0^{T}D\left(t\right)\sin \mathrm{n\omega tdt}(n=\mathrm{1,2},...\stackrel{-}{n})---(F1-4)$
$V_n=\frac{\omega }{\pi }{\int }_0^{T}U\left(t\right)\cos \mathrm{n\omega tdt}(n=\mathrm{0,1,2},...\stackrel{-}{n})---(F1-5)$
${\delta }_n=\frac{\omega }{\pi }{\int }_0^{T}U\left(t\right)\sin \mathrm{n\omega tdt}(n=\mathrm{1,2},...\stackrel{-}{n})---(F1-6)$
式中ω-曲柄角速度；T-抽汲周期。
实际工作中D(t)及U(t)是以曲线(或数值)形成给出的，所以傅里叶系数 可用近似的数值积分来确定。
以式(F1-1)和(F1-2)为边界条件，用分离变量法解方程(F1-0)，
可得抽油杆柱任意深度x断面的位移随时间的变化关系：
$U$$(x,t)=\frac{{\sigma }_o}{2E{A}_r}$$x+\frac{r_o}{2}$$\underset{n=1}{\overset{\stackrel{-}{n}}{\Sigma }}{[O}_n\left(x\right)\cos \mathrm{n\omega }$$+P_n\left(x\right)\sin \mathrm{n\omega t}$$]---$$(F1-7)$
根据虎克定律：
$F(x,t)=E{A}_r\frac{\partial U(x,t)}{\partial x}---(F1-8)$
则抽油杆柱任意深度x断面上的动载荷函数随时间的变化为：
$F(x,t)=E{A}_r[\frac{{\sigma }_o}{2E{A}_r}+\underset{n=1}{\overset{\stackrel{-}{n}}{\Sigma }}(\frac{\partial O_n\left(x\right)}{\partial x}\cos \mathrm{n\omega t}+\frac{\partial P_n\left(x\right)}{\partial x}\sin \mathrm{n\omega t})]---(F1-9)$
在t时间，x断面上的总载荷等于F(x，t)加x断面以下的抽油杆柱的重 量：
On(x)＝(Knchβnx+δnshβnx)sinanx+(μnshβnx+vnchβnx)cosanx   (F1-10)
Pn(x)＝(Knshβnx+δnchβnx)cosanx-(μnchβnx+vnshβnx)sinanx    (F1-11)
$a_n=\frac{\mathrm{n\omega }}{a\sqrt{2}}\sqrt{1+\sqrt{1+{\left(\frac{c}{\mathrm{n\omega }}\right)}^2}}---(F1-12)$
${\beta }_n=\frac{\mathrm{n\omega }}{a\sqrt{2}}\sqrt{-1+\sqrt{1+{\left(\frac{c}{\mathrm{n\omega }}\right)}^2}}---(F1-13)$
$K_n=\frac{{\sigma }_n{a}_n+{\tau }_n{\beta }_n}{E{A}_r(a_n^2+{\beta }_n^2)},{\mu }_n=\frac{{\sigma }_n{\beta }_n-{\tau }_n{a}_n}{E{A}_r(a_n^2+{\beta }_n^2)}---(F1-14)$
上述公式适用于单级抽油杆柱，对于多级抽油杆柱只需要做相应的扩充 就可得到类似的计算式。
对于多级杆来讲，根据两级杆相接点处的位移连续和载荷连续性，可以 得出付氏系数的递推公式。设在第i级杆内的位移和载荷为
$u_i(x,t)=a_{i0}x+b_{i0}+\underset{n=1}{\overset{\infty }{\Sigma }}[a_{\mathrm{in}}{e}^{(p_{\mathrm{in}}+q_{\mathrm{in}}j)x}+b_{\mathrm{in}}{e}^{-(p_{\mathrm{in}}+q_{\mathrm{in}}j)x}]e^{j{{\omega }_n}^t}---(F1-15)$
$F_i(x,t)={\left(\mathrm{EA}\right)}_i{a}_{i0}+\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\mathrm{EA}\right)}_i(P_{\mathrm{in}}+q_{\mathrm{in}}j)[a_{\mathrm{in}}{e}^{(p_{\mathrm{in}}+q_{\mathrm{in}}j)x}-b_{\mathrm{in}}{e}^{-{(p_{\mathrm{in}}+q_{\mathrm{in}}j)}^x}]e^{j{\omega }_{n}t}---(F1-16)$
则在第i级杆与第i+1级杆的交接点处x＝xi有
ui+1(xi，t)＝ui(xi，t) Fi+1(xi，t)＝Fi(xi，t)
(EA)i+1ai+1，0＝ai0·(EA)i
ai+1，0xi+bi+1，0＝ai0xi+bi0
(EA)i+1(Pi+1，n+qi+1，nj)ai+1，n＝(EA)i(Pin+qinj)ain
(EA)i+1(Pi+1，n+qi+1，nj)bi+1，n＝(EA)i(Pin+qinj)bin            (F1-17)
即
$a_{i+\mathrm{1,0}}=a_{i,0}\frac{{\left(\mathrm{EA}\right)}_i}{{\left(\mathrm{EA}\right)}_{i+1}}$
$b_{i+\mathrm{1,0}}=b_{i0}+a_{i0}(1-\frac{{\left(\mathrm{EA}\right)}_i}{{\left(\mathrm{EA}\right)}_{i+1}})x_i---(F1-18)$
$a_{i+1,n}=\frac{{\left(\mathrm{EA}\right)}_i(P_{\mathrm{in}}+q_{\mathrm{in}}j)}{{\left(\mathrm{EA}\right)}_{i+1}(P_{i+1,n}+q_{i+1,n}j)}{a}_{\mathrm{in}}$
$b_{i+1,n}=\frac{{\left(\mathrm{EA}\right)}_i(P_{\mathrm{in}}+q_{\mathrm{in}}j)}{{\left(\mathrm{EA}\right)}_{i+1}(P_{i+1,n}+q_{i+1,n}j)}{b}_{\mathrm{in}}$
其中a1n、b2n根据地面位移和载荷边界来确定，即
$a_{10}=\frac{F_0}{{\left(\mathrm{EA}\right)}_1},b_{10}=U_0---(F1-20)$
$a_{\ln }=\frac{1}{2}[U_n+\frac{F_n}{{\left(\mathrm{EA}\right)}_1}],b_{\ln }=\frac{1}{2}[U_n-\frac{F_n}{{\left(\mathrm{EA}\right)}_1}]---(F1-21)$
此外，还可采用杆柱波动方程的有限元差分解法。
根据出抽机井系统的工作过程所对应的力学行为，建立了其有限元差分 法模型，提出了有限差分的具体解法。数学模型波动方程中含有位移对时间t 的二阶偏导数和位移对位置x的二阶偏导数，一般来说，解这个波动方程将需 要两个初始条件和两个边界条件，然而在该预测模型中，仅仅需要两个边界 条件就够了，因为抽油杆的运动是具有周期性，况且，在实际中只需要一个 周期解，周期解与初始条件无关，无论把起始位置放到何处，其解均是相同 的。因此，也就不需要具体的初始条件。这样，实质上所研究的问题就是一 个包括偏微分方程和边界条件的一个边值问题。
偏微分方程难以求解析解，因此，对该数学模型进行数值求解。由于抽 油杆在井底的位移是未知待求的，用常规方法求解较为困难，可采用有限元 差分法，从井口开始，对每个单元输入i节点的载荷和位移，用单元方程求出 j节点的内力和位移，直至井底，求出泵处的内力和位移，得出泵功图。
差分法就是将微分算子(微分方程)转化为在节点的状态，这种方法称 为离散化，即要用差分方法求解某个问题，首先需要把自变量区域离散化。 具体讲，就是将位置空间(s)和时间空间(t)分割成若干个小区间，如果 分割成许多个小矩形，称此小矩形为网格。网格的顶点(xi，tj)称为节点。为 了将预测模型化为差分格式，将抽油杆柱离散成若干个单元，并以下标 i(i＝0，1，2，...，N)表示每个节点的位置，在时间方面取等步长Δt，以下标 j(j＝0，1，2，...M)表示，则ui，j表示抽油杆柱上第i个节点在时间tj时的位移。 根据差商概念可得：
对于第一级杆柱，除了两级杆界面点i以外各点差分格式如下：
${\left(\frac{\partial u}{\partial t}\right)}_{i,j}=\frac{u_{i,j+1}-u_{i,j}}{\mathrm{\Delta t}}+o\left(\mathrm{\Delta t}\right)---(F1-22)$
${\left(\frac{\partial u}{\partial t}\right)}_{i,j-1}=\frac{u_{i,j}-u_{i,j-1}}{\mathrm{\Delta t}}+o\left(\mathrm{\Delta t}\right)---(F1-23)$
${\left(\frac{{\partial }^{2}u}{\partial t^2}\right)}_{i,j}=\frac{{\left(\frac{\partial u}{\partial t}\right)}_{i,j}-{\left(\frac{\partial u}{\partial t}\right)}_{i,j-1}}{\mathrm{\Delta t}}=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j}}{\Delta t^2}+o\left(\Delta t^2\right)---(F1-24)$
${\left(\frac{\partial u}{\partial x}\right)}_{i,j}=\frac{u_{i+1,j}-u_{i,j}}{\mathrm{\Delta x}}+o\left(\mathrm{\Delta x}\right)---(F1-25)$
${\left(\frac{\partial u}{\partial x}\right)}_{i-1,j}=\frac{u_{i,j}-u_{i-1,j}}{\mathrm{\Delta x}}+o\left(\mathrm{\Delta x}\right)---(F1-26)$
${\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i,j}=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}+o\left(\Delta x^2\right)---(F1-27)$
将上述各式代入(F1-0)并略去余项得：
$\frac{u_{i.j+1}-2u_{i,j}+u_{i,j-1}}{\Delta t^2}={a_r}^2\left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}\right)-c_r\left(\frac{u_{i,j+1}-u_{i,j}}{\mathrm{\Delta t}}\right)+g_r---(F1-28)$
$u_{i.j+1}-2u_{i,j}+u_{i,j-1}=\frac{{a_r}^2\Delta t^2}{\Delta x^2}(u_{i+1,j}-2u_{i,j}+u_{i-1,j})-c_r\mathrm{\Delta t}(u_{i,j+1}-u_{i,j})+g_r\Delta t^2$
$u_{i.j+1}-2u_{i,j}+u_{i,j-1}={\lambda }_r(u_{i+1,j}-2u_{i,j}+u_{i-1,j})-c_r\mathrm{\Delta t}(u_{i,j+1}-u_{i,j})+g_r\Delta t^2$
式中： ${\lambda }_r=\frac{{a_r}^2\Delta t^2}{\Delta t^2}$
ui，j+1(1+crΔt)＝(2-2λr+crΔt)ui，j-ui，j-1+λr(ui+1，j+ui-1，j)+grΔt2
$u_{i.j+1}=\frac{(2-2{\lambda }_r+c_r\mathrm{\Delta t})}{(1+c_r\mathrm{\Delta t})}{u}_{i,j}-\frac{1}{(1+c_r\mathrm{\Delta t})}{u}_{i,j-1}+\frac{{\lambda }_r}{(1+c_r\mathrm{\Delta t})}(u_{i+1,j}+u_{i-1,j})+\frac{\Delta t^2}{(1+c_r\mathrm{\Delta t})}{g}_r$
$u_{i.j+1}=\frac{(2-2{\lambda }_r+c_r\mathrm{\Delta t})}{(1+c_r\mathrm{\Delta t})}{u}_{i,j}-\frac{1}{(1+c_r\mathrm{\Delta t})}{u}_{i,j-1}+\frac{{\lambda }_r}{(1+c_r\mathrm{\Delta t})}(u_{i+1,j}+u_{i-1,j})+\frac{\Delta t^2}{(1+c_r\mathrm{\Delta t})}{g}_r$
ui，j+1＝λr1ui，j-λr2ui，j-1+λr3(ui+1，j+ui-1，j)+λr4gr          (F1-29)
式中： ${\lambda }_{r1}=\frac{(2+2{\lambda }_r+c_r\mathrm{\Delta t})}{(1+c_r\mathrm{\Delta t})},$ ${\lambda }_{r2}=\frac{1}{(1+c_r\mathrm{\Delta t})},$ ${\lambda }_{r3}=\frac{{\lambda }_r}{(1+c_r\mathrm{\Delta t})},$ ${\lambda }_{r4}=\frac{\Delta t^2}{(1+c_r\mathrm{\Delta t})}$
第一级杆柱上除了i外任一点载荷的差分格式如下：
$P_{\left(r\right)i,j}=E_r{A}_r{\left(\frac{\partial u}{\partial x}\right)}_{i,j}=E_r{A}_r\frac{u_{i+1,j}-u_{i-1,j}}{2\mathrm{\Delta x}}---(F1-30)$
同理第二级常规杆(加重杆)的差分格式除了i点以外各点的差分格式如下：
ui，j+1＝λj1ui，j-λj2ui，j-1+λj3(ui+1，j+ui-1，j)+λrjgj     (F1-31)
式中： ${\lambda }_{j1}=\frac{2-2{\lambda }_j+c_j\mathrm{\Delta t}}{(1+c_j\mathrm{\Delta t})},$ ${\lambda }_{j2}=\frac{1}{(1+c_j\mathrm{\Delta t})},$
${\lambda }_{j3}=\frac{{\lambda }_j}{(1+c_j\mathrm{\Delta t})},$ ${\lambda }_{j4}=\frac{\Delta t^2}{(1+c_j\mathrm{\Delta t})};$ ${\lambda }_j=\frac{{a_j}^2\Delta t^2}{\Delta x^2}$
第二级杆柱上除了i点外各点载荷的差分格式如下：
$P_{\left(j\right)i,j}=E_j{A}_j{\left(\frac{\partial u}{\partial x}\right)}_{i,j}=E_j{A}_j\frac{u_{i+1,j}-u_{i-1,j}}{2\mathrm{\Delta x}}---(F1-32)$
对于钢丝绳和加重杆的连接部分，由于钢丝绳和加重杆的连接长度很小， 所以对整个杆体的动力学特性影响很小，这里采用等效值的方法处理：
$E^{\prime }=\frac{L^2}{(\frac{L_r}{E_r{A}_r}+\frac{L_j}{E_j{A}_j})(L_r{A}_r+L_j{A}_j)}$
${\rho }^{\text{'}}=\frac{{\rho }_r{A}_r+{\rho }_j{A}_j}{(A_r+A_j)};$ $A^{\text{'}}=\frac{L_r{A}_r+L_j{A}_j}{L}$
差分格式同(F1-34)
另外，根据地面示功图计算井下示功图时，必须首先确定阻尼系数。抽 油杆柱系统的阻尼力包括粘阻尼力和非粘滞阻尼力。粘滞阻尼力有抽油杆、 接筛与液体之间的粘滞摩擦力、泵阀和阀座内孔的流体压力损失等。
钢丝绳杆既没有接箍，在绳体的表面还有防腐涂层。因此，在相同条件 下，钢丝绳杆在流体中的阻尼远比常规钢杆的小。计算时，可用由抽油杆柱 在一个循环中由粘滞阻尼引起的摩擦功来确定的阻尼系数：
$c_r=\frac{2\mathrm{\pi \eta }}{{\rho }_r{A}_r}f\left(m\right)---(F1-33)$
其中： $f\left(m\right)=\frac{-2\ln m+m^2-1}{2(\ln m+m^2\ln m-1+m^2)}+\frac{1}{\ln m}$
考虑到钢丝绳的表面结构特性，取无量纲修正系数ψ
则 $c_r=\psi \frac{2\mathrm{\pi \eta }}{{\rho }_r{A}_r}f\left(m\right)---(F1-34)$
式中：
dr-抽油杆直径；dt-油管内径； $m=\frac{d_t}{d_r};$ η-液体粘度，Pa·s；ρr-抽油 杆的密度，Kg/m3；Ar-抽油杆的截面积，m2；
最终以井下泵功图和油井的气油比、含水、油粘度等物理参数的变化为 依据，对井下泵工作状况进行诊断和各项指标的量化，计算出泵的有效冲程、 充满系数、气影响程度，计算泵的有效排量，进而折算求出井口有效排量。
实际应用中，通过对计算产液量值的标定，可以获得较高的计量准确度。 抽油机井泵功图识别得到的是泵工作压力下的地下液体体积流量，折算求出 井口有效排量主要是根据质量连续原理，应用原油压缩系数的，一个油田原 油的压缩系数是比较固定的，一般都约等于1，不会偏离1很大的距离，所以 一般地面体积和地下体积相差不多，在没有原油压缩系数的情况下可以将这 个系数涵盖于地面流量标定系数中。
2.自喷井
如图6所示，为自喷井时计算产液量的流程图。其中，应用节流压降计算 模型应用节流压差计算液体的流过量。该数学模型中考虑了气体在不同状况 下所占体积。本模型中，根据的工况数据指：瞬时产液量、累积产液量等； 根据的基础数据指：嘴前嘴后压力P1，P2等动态参数，油嘴直径d、生产气液比Rs 等静态参数。
其中，自喷井产液量计算主要根据动态参数：嘴前压力P1(油压)、嘴后 压力P2(回压)；静态参数：油嘴直径d、生产气液比Rs，结合多相流油嘴节 流模型计算出自喷井体积流量，然后用流量标定系数进行修正，得到自喷井 井口标准流量，并计算出油井的流压与产能。
其理论模型示意函数如下：
Q自喷＝kf(d，P1，P2，Rs)                                         (F2-0-1)
Q自喷----自喷井的产液量，m3；d----油嘴直径，mm；P1----嘴前压力(油 压)，MPa；P2----嘴后压力(回压)，MPa；Rs----生产气液比；k----流量 标定系数。
应用贝努利方程和流动连续性原理，可推导出以下的理论流量公式，推 导过程从略：
$q_v=\frac{1}{\sqrt{1-{\beta }^4}}\frac{\pi }{4}{d}^2\sqrt{\frac{2\mathrm{\Delta p}}{{\rho }_1}}---(F2-0-2)$
上式适用于不含气体的纯液体计算，如果油井含有天然气，需要进行天 然气校正，又由于校正系数c的定义是：c＝实际流量/理论流量，就是前面的流 量标定系数，最后可得出油嘴节流适用的流量公式：
$q_v=\frac{c·ϵ}{\sqrt{1-{\beta }^4}}\frac{\pi }{4}{d}^2\sqrt{\frac{2\mathrm{\Delta P}}{{\rho }_1}}(F2-0-3)$
ε的计算公式如下：
$ϵ=1-(0.649+0.696{\beta }^4)·\frac{\mathrm{\Delta P}}{k·p_1}---(F2-0-4)$
质量流量qm＝qv·ρ1
以上各式中：
ε——被测介质的可膨胀性系数，对于液体ε＝1；对气体、蒸气等可压缩 流体ε＜1；qv——流体的体积流量，[m3/s]；(工况下总流体的体积流量)； qm——总流体的质量流量，[kg/S]；d——工作状况下节流件的等效开孔直径， [m](对于孔板是孔径，对于文丘得管是喉径)；ΔP-节流差压，即油压-回 压，ΔP＝P1-P2；[Pa]；ρ1——工作状况下，节流前上游处流体的密度，[kg/m3]；
根据气油比等参数计算出，这里从略；C——流出系数，无量纲；β——直径 比，无量纲，与节流嘴的形状有关，对于锥形的节流嘴，β＝d/D，即小直径/ 大直径，对于一般通用的油嘴为等直径的孔板，β＝1；k-天然气的等熵指数， 通常取1.229；P1-工况下节流件(内锥)上游取压孔处可压缩流体的绝对静 压，既油井的油压Pa。
另外，自喷井部分的多相流计算部分如下所述：
其中，F2.1多相流计算的基本参数如下：
(1)多相流压力梯度方程
多相管流的压力梯度包括：因举高液体而克服重力所需的压力势能、流 体因加速而增加的动能和流体沿管路的摩阻损失，其数学表达式如下：
$\frac{\mathrm{dp}}{\mathrm{dh}}={\rho }_{m}g\sin \theta +{\rho }_m{v}_m\frac{{\mathrm{dv}}_m}{\mathrm{dh}}+f_m\frac{{\rho }_m}{d}\frac{v_m^2}{2}---(F2-1)$
式中ρm为多相混合物的密度；vm为多相混合物的流速；fm为多相混合物 流动时的摩擦阻力系数；d为管径；p为压力；h为深度；g为重力加速度；θ为 井斜角的余角。
(2)摩阻计算
此问题是已知管道内的流量Q、管道内径d、管道长度L、水的密度ρ、水 的动力粘度μ、管道内壁的绝对粗糙度Δ，求此段管道中的压力损耗。此问 题可按下面步骤计算。
根据管道内径由下式计算管道的内截面积A：
$A=\frac{\pi }{4}{d}^2---(F2-2)$
根据管道内的流量和管道的内截面积计算管道内水流的平均流速v：
v＝Q/A                                                         (F2-3)
根据管道内水流的平均流速、水的密度、水的粘度和管道内径计算管道 内水流的雷诺数Re：
$\mathrm{Re}=\frac{\mathrm{vd\rho }}{\mu }---(F2-4)$
根据管道的绝对粗糙度和管道内径计算管道内壁的相对粗糙度ε：
$ϵ=\frac{2\Delta }{d}---(F2-5)$
根据管道内壁的相对粗糙度和雷诺数从表F2-2所列公式中选取一个合适 用公式计算摩阻系数λ。
根据摩阻系数、管道长度、管道内径、管道内水流速度和水的密度由下 式计算压力损耗ΔP。
$\mathrm{\Delta P}=\lambda \frac{L}{d}\rho \frac{v^2}{2}---(F2-6)$
当途中有分流时，沿程流量改变，应分段进行计算。
摩阻系数计算公式见表1。
表1

(3)按深度增量迭代的步骤
根据多相垂直管流的压力梯度就可计算出沿程压力分布。由于多相垂直 管流中每相流体的物理参数及混合物密度和流速都随压力和温度而变，沿程 压力梯度并不是常数。因此，多相管流需要分段计算，并要预先求得相应段 的流体性质参数。然而，这些参数又是压力和温度的函数，压力却又是计算 中需要求的未知数。所以多相管流需要采用迭代法进行计算。下面是通用的 计算步骤：
a.已知任一点的压力P0作为起点，任选一个合适的深度间隔Δh。
b.估计一个对应于计算间隔的压力增量ΔP。
c.计算本段的平均温度 和平均压力 ，并确定该 知 下的全部流体性 质参数：溶解油气比Rs、原油体积系数B0和粘度μ0、气体密度ρg和粘度μg、 混合物粘度μm、液体的表面张力σ。
d.先计算该段的压力梯度 ，然后计算对应于Δh的压力增量 $\mathrm{\Delta p}=\mathrm{\Delta h}{\left(\frac{\mathrm{dp}}{\mathrm{dh}}\right)}_i.$
e.比较估计和计算的压力增量，若二者之差不在允许范围之内，则以计 算值作为新的估计值，重复第(2)至(5)步，直到两者之差在允许范围之 内为止。
f.计算该段下端对应的深度Li和压力pi：Li＝iΔh； $P_i=P_0+\underset{i=1}{\overset{n}{\Sigma }}\Delta p_i;$ f=1，2，…，n。
g.以Li处的压力pi为起点压力重复第(2)至(7)步，计算下一段的深 度Li+1和压力Pi+1，直到各段累加深度等于或大于管长L(Ln≥L)时为止。
(4)环形流道的当量直径
流体在环形流道中的摩阻可近似地用圆管的公式。此时，按照圆管与环 形流道的水力斗径相等，把环形流道化为摩阻相当的圆管井行计算。即
$R=\frac{d_e}{4}---(F2-7)$
式中R为水力半径；de为当量直径。
环形流道的面积为：
$A=\frac{\pi }{4}(d_c^2-d_t^2)---(F2-8)$
式中dc为套管的内径；dt为套管的外径。
环形流道的湿周为：
x＝π(dc+dt)                                                 (F2-9)
因此，根据水力半径的定义，可求得环形流道的水力半径：
$R=\frac{A}{\chi }=\frac{(d_c-d_t)}{4}---(F2-10)$
环形流道的当量直径为：
de＝4R＝dc-dt                                                (F2-11)
F2.2 Hagedorn-Brown模型
Hagedorn-Brown(1965)基于所假设的压力梯度模型，根据大量的现场 实验数据反算持液率，提出了应用于各种流型下的两相垂直上升管流压降关 系式，此压降关系式不需要判别流型适用于产水气井的流动条 件，Hagedorn-Brown法适用于低气液比的高产排液井。
(1)压力梯度方程式
$\frac{\mathrm{\Delta p}}{\mathrm{\Delta h}}=g{\rho }_m\sin \theta +\frac{f_m{\rho }_n^2{v}_m^2}{2d{\rho }_m}+{\rho }_m\frac{\Delta \left(\frac{v_m^2}{2}\right)}{\mathrm{\Delta h}}---(F2-12)$
(2)计算方法
①判断流型
Hagedorn-Brown模型不需判别流型，只是当Nvg＜L1时使用Griffith相关 式计算混合物密度ρm和fm，具体细节见Orkiszewski方法。其中
$L_1=1.071-0.7277v_m^2/d\ge 0.13$
②计算混合物的平均密度
a.由图10确定无因次液相粘度准数Nl；
b.计算滞留量的相关函数U：
$U=N_{\mathrm{vl}}{(p/0.101)}^{0.1}{N}_l{N}_{\mathrm{vg}}^{-0.575}/N_d---(F2-13)$
由图11求出φ/ψ；
c.由图12求出ψ；
d.计算持液率：
φ′＝φ′l/ψ)ψ                                           (F2-14)
e.计算混合物的平均密度ρm：
(a)由φ′计算ρm(φ)
ρm(φ′)＝ρlφ′+ρg(1-φ′)                              (F2-15)
(b)比较ρm(φ′)和ρm，采用其中较大的值。
③确定摩阻系数fm
a.计算两相雷诺数NREm：
$N_{\mathrm{REm}}=1000{\rho }_n{v}_{m}d/\left({\mu }_l^{{\phi }^{\text{'}}}{\mu }_g^{1-{\phi }^{\text{'}}}\right)---(F2-16)$
b.根据ε/d和NREm，由公式确定摩阻系数fm。
④计算Δ(Vm2)
F2.3Orkiszewski方法
Orkiszewski法适用于高气油比的中低产井。
(1)压力梯度方程式
$\frac{\mathrm{\Delta p}}{\mathrm{\Delta h}}=\frac{g{\rho }_m+{\tau }_f}{1-{\rho }_m{v}_m{v}_{\mathrm{sg}}/p}---(F2-17)$
(2)计算方法
①判断流型
    界限    qsg/qm＜L1    qsg/qm＞L1，       Nvg＜L2    L3＞Nvg＞       L2    Nvg＞L3     流型      气泡流      段塞流      过渡流      环雾流
式中：
$L_1=1.071-0.7277v_m^2/d\ge 0.13$
L2＝50+36Nvgqsl/qsg
L3＝75+84(Nvgqsl/qsg)0.75
②根据流动方式，确定具体的平均密度ρm和摩阻损耗梯度τf。
a.气泡流
(a)气体所占的空间分数Hg：
$H_g=\{1+\frac{q_m}{v_{s}A}-{[{(1+\frac{q_m}{v_{s}A})}^2-\frac{4q_{\mathrm{sg}}}{v_{s}A}]}^{0.5}\}/2---(F2-18)$
式中vs取0.244m/s是一个较好的近似值。
(b)流体平均密度ρm：
ρm＝(1-Hg)ρl+Hgρg                                            (F2-19)
(c)摩阻梯度τf：
${\tau }_f=f{\rho }_l{v}_{\mathrm{sl}}^2/\left[2d(1-H_g)\right]---(F2-20)$
式中f根据ε/d和NRE由公式计算，NRE＝1000ρlvsld/[μl(1-Hg)]。
b.段塞流
(a)流体平均密度ρm：
ρm＝(wm+ρlvsA)/(qm+vsA)+δρl                                 (F2-21)
式中滑脱速度vs是气泡流雷诺数NREb和雷诺数NRE的函数：
NREb＝1000ρlvsd/μl
NRE＝1000ρlvmd/μl
当NREb≤3000时：
vs＝(0.546+8.74×10-6NRE)(gd)0.5                                (F2-22)
当NREb≥8000时：
vs＝(0.35+8.74×10-6NRE)(gd)0.5                                 (F2-23)
当8000＞NREb3000时：
vs0＝(0.251+8.74×10-6NRE)(gd)0.5                               (F2-24)
$v_s=\{v_{s0}+{[v_{s0}^2+0.1202{\mu }_l/\left({\rho }_l{d}^{0.5}\right)]}^{0.5}\}/2---(F2-25)$
液相分布系数δ与及连续液相有关：
a)油为连续相，并且vm＞3.048m/s时，
δ＝0.00537lg(μ1+1)/d1.371+0.569lg(d)+0.455-lg(vm/0.3048)
[0.001574lg(μ1+1)/d1.571+063lg(d)+0.722)]                      (F2-26)
b)水为连续相，并且vm＞3.048m/s时，
δ＝0.01741g(μ1)/d0.799-0.888lg(d)-0.162lg(vm)-1.2508          (F2-27)
c)油为连续相，并且vm＜3.048m/s时，
δ＝0.024lg(μ1+1)/d1.415+0.1131g(d)+0.1671g(vm)-1.1395          (F2-28)
d)水为连续相，并且vm＜3.048m/s时，
δ＝0.00252lg(μ1)/d1.35-0.428lg(d)+0.2321g(vm)-0.782           (F2-29)
δ还要受以下条件的限制：
当vm＜3.048m/s时，δ≥-0.2133vm；
当vm＞3.048m/s时，δ≥-vsAp(1-ρm/ρl)/(qm+vsA)。
(b)摩阻梯度τf：
${\tau }_f=\frac{f{\rho }_l{v}_m^2}{2d}(\frac{q_{\mathrm{sl}}+v_{s}A}{q_m+v_{s}A}+\delta )---(F2-30)$
式中f根据ε/d和NRE由公式计算，NRE＝1000ρlvmd/μl。
c.过渡流
过渡流流型的平均密度ρm和摩阻损耗梯度τf，均按Nvg进行段塞流和环雾 流的线性加权平均：
ρm＝[(Lm-Nvg)ρm段+(Nvg-Ls)ρm雾]/(Lm-Ls)                      (F2-31)
τf＝[(Lm-Nvg)τf段+(Nvg-Ls)τf雾]/(Lm-Ls)                      (F2-32)
为了更准确地预测摩阻损耗梯度，环雾流中的qsg按下式计算：
$q_{\mathrm{sg}}=A{L}_m{\left(\frac{{\rho }_l}{g{\sigma }_l}\right)}^{0.25}---(F2-33)$
d.环雾流
(a)气体所占的空间分数Hg：
Hg＝qsg/qm                                                      (F2-34)
(b)流体的平均密度ρm：
ρm＝(1-Hg)ρl+Hgρg                                            (F2-35)
(c)摩阻梯度τf：
${\tau }_f=f{\rho }_g{v}_{\mathrm{sg}}^2/\left(2d\right)---(F2-36)$
式中f根据ε/d和NRE由公式计算，NRE＝1000ρgvsgd/μg。
由于环雾流中有液膜形成，相对粗糙度需重新计算，但它受条件0.001＜ ε/d＜0.5的限制。ε/d的计算方法为：
Nw＜0.005时：
$ϵ/d=0.01543{\sigma }_l/\left({\rho }_g{v}_{\mathrm{sg}}^{2}d\right)---(F2-37)$
Nw≥0.005时：
$ϵ/d=0.08078{\sigma }_l{N}_w^{0.302}/\left({\rho }_g{v}_{\mathrm{sg}}^{2}d\right)---(F2-38)$
式中Nw的定义为：
Nw＝4.865×10-6(vsgμl/σl)2(ρg/ρl)                           (F2-39)
F2.4  Beggs-Brill方法
Beggs-Brill方法适用于各种角度的管线，井中油管和环空管，斜井时建 议使用。
(1)压力梯度方程式
$\frac{\mathrm{\Delta p}}{\mathrm{\Delta h}}=\frac{g{\rho }_m\sin \theta +f_m{\rho }_n{v}_m^2/\left(2d\right)}{1-{\rho }_m{v}_m{v}_{\mathrm{sg}}/p}---(F2-40)$
(2)计算方法
①判段流型
                         界限    流型      λ＜0.01及NFR＜L1或λ≥0.01及NFR＜L2    分异流         λ≥0.01及L2＜NFR≤L3    过渡流     0.01≤λ＜0.4及L3＜NFR≤L1或λ≥0.4及L3＜NFR≤L4    间隔流      λ＜0.4及NFR≥L1或λ≥0.4及NFR＞L4    分散流
式中：
L1＝316λ0.302；L2＝0.0009252λ-2.4684；L3＝0.1λ-1.4516；L4＝0.5λ-6.738
②计算持液率Hl及两相密度ρM
a.计算水平持液率Hl(0)：
$H_l\left(0\right)=G_1{\lambda }^{C_2}/N_{\mathrm{FR}}^{C_3}---(F2-41)$
式中回归系数C1，C2和C3的取值见表2。
表2水平持液率公式回归系数
     流型      C1      C2       C3     分异流     0.98     0.4846     0.0868     间隔流     0.84     0.5 351     0.0173     分散流     1.065     0.5824     0.0609
b.计算斜度校正系数
$C=(1-\lambda )\ln \left(C_4{\lambda }^{C_5}{N}_{\mathrm{vl}}^{C_6}{N}_{\mathrm{FR}}^{C_7}\right)---(F2-42)$
式中回归系数C4，C5，C6和C7的取值见表3。
表3斜度校正系数公式回归系数
       流型       C4       C5       C6       C7      分异流，上坡     0.011     -3.768     3.539     -1.614      间隔流，上坡     2.96     0.305     -0.4473     0.0978      分散流，上坡          不必校正C＝0      各种流型，下坡     4.70     -0.3692     0.1244     -0.5056
c.计算持液率校正系数ψ：
ψ＝1+C[sin(1.8θ)-0.333sin3(1.8θ)]                            (F2-43)
对于垂向井：ψ＝1+0.3C
d.计算持液率Hl(θ)及两相密度ρm：
Hl(θ)＝Hl(0)ψ
ρm＝ρlHl+ρg(1-Hl)                                            (F2-44)
③计算两相摩阻系数fm：
a.计算摩阻系数比(fm/fn)：
(fm/fn)＝es                                                     (F2-45)
式中：
S＝ln(y)/{-0.0523+3.182ln(y)-0.8725[ln(y)]2+0.01853[ln(y)]4}
y＝λ/[Hl(θ)]2
当1＜y＜1.2时：s＝ln(2.2y-1.2)
b.计算无滑脱的摩阻系数fn：
fn＝{21g[NREn/(4.5223lgNREn-3.8215)]}-2                         (F2-46)
或 $f_n=0.0056+0.5/N_{\mathrm{REn}}^{0.32}---(F2-47)$
c.计算两相摩阻系数fm：
F2.5 Hasan-Kabir(哈森)方法
哈森方法将流动分为四种流型：泡流、段塞流、涡流和环流。对于不同 的流型，采用不同的方法计算流体参数和由于摩阻产生的压力梯度。
气泡在垂直系统中最终上升速度的表达式可由气泡在液体中运动时的浮 力和阻力的平衡导出。
$v_{\infty }=1.53{\left[\frac{\mathrm{g\sigma }({\rho }_L-{\rho }_g)}{{\rho }_L^2}\right]}^{0.25}---(F2-48)$
式中v∞为气泡的最终上升速度；g为重力加速度；σ为液体的表面张力；ρL 为液体的密度；ρg为气体的密度。
Taylor泡的上升速度为：
$v_{\infty T}=0.35\sqrt{\frac{\mathrm{gd}({\rho }_L-{\rho }_g)}{{\rho }_L}}---(F2-49)$
式中v∞T为气泡上升的Taylor速度；d为流道的直径。
(1)泡流
当流动的流体满足式(F2-50)或(F2-51)时，则流体流动的流型为泡流。
此时，摩阻产生的压力梯度由式(F2-54)计算。
vgs＜0.429vLs+0.357v∞                                          (F2-50)
 式中vgs为气体的表观速度，即气体流量除以整个流道的截面积；vLs为液 体的表观流速。
fg＜0.52且 $V_m^{1.12}>4.68d^{0.48}{\left[\frac{g({\rho }_L-{\rho }_g)}{\sigma }\right]}^{0.5}{\left(\frac{\sigma }{{\rho }_L}\right)}^{0.6}{\left(\frac{{\rho }_m}{{\mu }_L}\right)}^{0.08}---(F2-51)$
式中fg为空隙因数；vm为气液混合物的流速；μL为液体的粘度。而空隙 因数由下列式子计算：
$f_g=\frac{v_{\mathrm{gs}}}{C_0{v}_m+v_{\infty }},$ $C_0=1.2+0.371\left(\frac{d_t}{d_c}\right)---(F2-52)$
式中C0为无量纲参数；dt为油管的外径；dc为套管的内径。
气液混合物的密度由下式计算：
ρm＝(1-fg)ρL+fgρg                                            (F2-53)
式中ρm为气液混合物的密度。
如流体的流型为泡流，则由摩阻产生的压力梯度由下式计算：
${\left(\frac{\mathrm{dp}}{\mathrm{dz}}\right)}_F=\frac{2f_m{v}_m^2{\rho }_m}{d}---(F2-54)$
式中的摩阻系数fm是雷诺数的函数：
$f_m\propto \left(\frac{{\mathrm{dv}}_m{\rho }_L}{{\mu }_L}\right)---(F2-55)$
(2)段塞流
当流道中的气液两相流动满足下列条件时，则流动的流型为段塞流。
vgs＞0.429vLs+0.357v∞                                          (F2-56)
且当 ${\rho }_L{v}_{\mathrm{Ls}}^2>50$ 时，
${\rho }_g{v}_{\mathrm{gs}}^2<17.11{\mathrm{og}}_{10}\left({\rho }_L{v}_{\mathrm{Ls}}^2\right)-23.2;---(F2-57)$
当 ${\rho }_L{v}_{\mathrm{Ls}}^2<50$ 时，
${\rho }_g{v}_{\mathrm{gs}}^2<0.00673{\left({\rho }_L{v}_{\mathrm{Ls}}^2\right)}^{1.7}.---(F2-58)$
当流体的流型为段塞流时，其空隙因数由下式计算：
$f_g=\frac{v_{\mathrm{gs}}}{C_1{v}_m+v_{\infty T}},$ $C_1=1.18+0.90\left(\frac{d_t}{d_c}\right)---(F2-59)$
当流体的流型为段塞流时，其密度由下式计算：
ρm＝(1-fg)ρL+fgρg                                            (F2-60)
当流体的流型为段塞流时，其中由于摩阻产生的压力梯度由下式计算：
${\left(\frac{\mathrm{dp}}{\mathrm{dz}}\right)}_F=\frac{2f_m{v}_m^2{\rho }_L(1-f_g)}{d}---(F2-61)$
式中的摩阻系数是雷诺数的函数：
$f_m\propto \left(\frac{{\mathrm{dv}}_m-{\rho }_L}{{\mu }_L}\right)---(F2-62)$
(3)涡流
当流道中的气液两相流动满足下列条件时，则流动的流型为涡流。
$v_{\mathrm{gs}}<3.1{\left[\frac{\mathrm{\sigma g}({\rho }_L-{\rho }_g)}{{\rho }_g^2}\right]}^{0.25}---(F2-63)$
且当 ${\rho }_L{v}_{\mathrm{Ls}}^2>50$ 时，
${\rho }_g{v}_{\mathrm{gs}}^2>17.11{\mathrm{og}}_{10}\left({\rho }_L{v}_{\mathrm{Ls}}^2\right)-23.2;---(F2-64)$
当 ${\rho }_L{v}_{\mathrm{Ls}}^2<50$ 时，
${\rho }_g{v}_{\mathrm{gs}}^2>0.00673{\left({\rho }_L{v}_{\mathrm{Ls}}^2\right)}^{1.7}.---(F2-65)$
当流动的流型为涡流时，其空隙因数由下式计算：
$f_g=\frac{v_{\mathrm{gs}}}{C_1{v}_m+v_{\infty T}},$ $C_1=1.15+0.90\left(\frac{d_t}{d_c}\right)---(F2-66)$
式中C1为无因次量。
当流动的流型为涡流时，其密度由下式计算：
ρm＝(1-fg)ρL+fgρg                                            (F2-67)
当流动的流型为涡流时，其中由摩阻产生的压力梯度由下式计算：
${\left(\frac{\mathrm{dp}}{\mathrm{dz}}\right)}_F=\frac{2f_m{v}_m^2{\rho }_L(1-f_g)}{d}---(F2-68)$
式中的摩阻系数是雷诺数的函数：
$f_m\propto \left(\frac{{\mathrm{dv}}_m{\rho }_L}{{\mu }_L}\right)---(F2-69)$
(4)环流
当流道中的气液两相流满足下式时，则其流型为环流。
$v_{\mathrm{gs}}>3.1{\left[\frac{\mathrm{\sigma g}({\rho }_L-{\rho }_g)}{{\rho }_g^2}\right]}^{0.25}---(F2-70)$
当流动的流型为环流时，其空隙因数由下式计算：
fg＝(1+X0.8)-0.378                                              (F2-71)
式中X为无因次量，它由下式计算：
$X={\left(\frac{{\rho }_g}{{\rho }_L}\right)}^{0.5}{\left(\frac{1-x}{x}\right)}^{0.9}{\left(\frac{{\mu }_L}{{\mu }_g}\right)}^{0.1}---(F2-70)$
式中x为气体的质量比例。
当流动的流型为环流时，其密度由下式计算：
${\rho }_c=\frac{v_{\mathrm{gs}}{\rho }_g+{\mathrm{Ev}}_{\mathrm{Ls}}{\rho }_L}{v_{\mathrm{gs}}+{\mathrm{Ev}}_{\mathrm{gs}}}---(F2-73)$
式中E为核心气体所携带的液体的比例。
当流动的流型为环流时，其中由于摩阻产生的压力梯度由下式计算：

式中NReg为雷诺数。
本油嘴节流计算模型同时还适用于使用地面油嘴的抽油机井、电潜泵井、 螺杆泵井的只知道地面节流前后压力的产液量计算，可以作为独立的方法应 用，也可以作为辅助的方法与主体的方法进行对比，考虑到模型的先进性， 本发明中抽油机井、电潜泵井、螺杆泵井的主要计算方法为本实施例的算法。
3.螺杆泵井
如图7所示，为螺杆泵井产液量计算流程图。根据关系式来计算出质量流 量。该关系式为：有功功率＝质量流量×排出压力。
所述工况数据为瞬时产液量、累积产液量，井口部分的扭矩、转速、功 率，泵端的扭矩、转速、功率，有效功率，系统损耗、系统效率、泵排出压 力，电量、电费等，还有丰富的宏观控制图、诊断结论等诸多指标。所述油 井基础数据为转速S、三相电流I1，I2，I3电压U、功率因素cosφ、扭矩M、载 荷P等动态参数；生产气液比Rs等静态参数根据。
如图7所示，主要根据动态参数：转速s、三相电流I1，I2，I3、电压U、功 率因素cosφ、扭矩、载荷P；静态参数：生产气液比Rs；利用力学计算数学 模型和功耗计算数学模型拟合，计算出经过流量标定系数k修正，计算出螺 杆泵井的地面标准状况下的产液量。其理论模型示意函数如下：
Q螺杆泵＝kf(S，I1，I2，I3，U，cosφ，M，P，Rs)                  (F3-0-1)
Q螺杆泵----螺杆泵的产液量，m3；s----转速，转/分；----扭矩，N·m；
----载荷，kN；Rs----生产气液比；I1，I2，I3----三相电流，A；U----电压， V；cosφ----功率因数；k----流量标定系数，小数。
如图8A和图8B所示，梁的动力学方程理论示意图。其中，螺杆泵井根据 梁的动力学方程理论如下：
假设梁不受外力作用，则梁的势能泛函为：
$\Pi =l{\int }_0^1\frac{1}{2}{ϵ}^T\mathrm{Dϵd\zeta }---(F3-0-2)$
动能泛函为
$\Gamma =l{\int }_0^1[\frac{1}{2}\mathrm{\rho A}({\stackrel{.}{u}}^2\stackrel{.}{+}{\stackrel{.}{v}}^2+{\stackrel{.}{w}}^2)+\frac{1}{2}\rho I_{\zeta }{{\stackrel{·}{\theta }}_{\zeta }}^2]\mathrm{d\zeta }=l{\int }_0^1\frac{1}{2}{\stackrel{·}{u}}^T\rho \stackrel{·}{u}\mathrm{d\zeta }---(F3-0-3)$
其中 $\rho =\left[\begin{array}{cccc}\mathrm{\rho A}& 0& 0& 0\\ 0& \mathrm{\rho A}& 0& 0\\ 0& 0& \mathrm{\rho A}& 0\\ 0& 0& 0& \rho I_{\zeta }\end{array}\right],$ ρ为梁的密度，Iζ为绕ζ轴的转动惯量。
根据Hamilton变分原理：
$\delta {\int }_{t_1}^{t_2}{\int }_0^1(\Pi -\Gamma )\mathrm{dldt}=0---(F3-0-4)$
可得梁的运动微分方程如下：
$\left.\begin{array}{c}\mathrm{EA}\frac{{\partial }^{2}u}{l^2\partial {\zeta }^2}-\mathrm{\rho A}\frac{{\partial }^{2}u}{\partial t^2}=0\\ {\mathrm{EI}}_{\eta }\frac{{\partial }^{4}v}{l^4\partial {\zeta }^4}+\mathrm{\rho A}\frac{{\partial }^{2}v}{\partial t^2}=0\\ {\mathrm{EI}}_{\xi }\frac{{\partial }^{4}w}{l^4\partial {\zeta }^4}+\mathrm{\rho A}\frac{{\partial }^{2}w}{\partial t^2}=0\\ {\mathrm{GJ}}_{\rho }\frac{{\partial }^2{\theta }_{\zeta }}{l^2\partial {\zeta }^2}-{\mathrm{\rho I}}_{\zeta }\frac{{\partial }^2{\theta }_{\zeta }}{\partial t^2}=0\end{array}\right\}---(F3-0-5)$
本发明中螺杆泵的适时参数主要计算就是通过求解以上方程组得到的， 以上波动方程的细节如下述的螺杆泵井杆柱动态模型：
将螺杆泵井杆柱作为一个2节点结构单元梁，引入位移的插值形式u＝Nae
其中ae为单元节点位移矢量，N为插值函数矩阵
ae＝(u1，v1，w1，θζ1，θξ1，θη1，u2，v2，w2，θζ2，θξ2，θη2)T    (F3-1)
$N=\left[\begin{array}{cccccccccccc}{N}_1& 0& 0& 0& 0& 0& N_2& 0& 0& 0& 0& 0\\ 0& H_1^{\left(0\right)}& 0& 0& 0& H_1^{\left(1\right)}& 0& H_2^{\left(0\right)}& 0& 0& 0& H_2^{\left(1\right)}\\ 0& 0& H_1^{\left(0\right)}& 0& -H_1^{\left(1\right)}& 0& 0& 0& H_2^{\left(0\right)}& 0& -H_2^{\left(1\right)}& 0\\ 0& 0& 0& N_1& 0& 0& 0& 0& 0& N_2& 0& 0\end{array}\right]---(F3-2)$
其中Ni为Lagrange插值函数，Hi(j)(i，j＝1，2)为Hermite插值函数
$\left.\begin{array}{c}N1=F3-\zeta ,N2=\zeta \\ H_1^{\left(0\right)}=1-3{\zeta }^2+2{\zeta }^3;H_1^{\left(1\right)}=(\zeta -2{\zeta }^2+{\zeta }^3)l\\ H_2^{\left(0\right)}=3{\zeta }^2-2{\zeta }^3,\\ H_2^{\left(1\right)}=({\zeta }^3-{\zeta }^2)l\end{array}\right\}---(F3-3)$
(F3-2)式中对轴向位移u和扭角θζ用Lagrange插值函数，对挠度v，w 采用Hermite插值函数，这样得到的节点处转角θη，θξ即分别为v，w的一阶 导数。将位移插值形式代入(F3-0-2)式的动能泛函与(F3-0-3)式的势能泛函 中，再由拉格朗日方程：
$\frac{d}{\mathrm{dt}}\left[\frac{\partial \Gamma }{\partial \stackrel{·}{u}}\right]-\frac{\partial \Gamma }{\partial u}+\frac{\partial \Pi }{\partial u}=0---(F3-4)$
可得梁单元的动力学有限元方程：
Mee+Keae＝0                                                   (F3-5)
式中Me为单元质量矩阵，Ke为单元刚度矩阵，并与预轴力有关，PTe为单元 等效温度载荷矢量。并有
$M^e=l{\int }_0^1{N}^T\mathrm{\rho Nd\zeta }---(F3-6)$
$K^e=l{\int }_0^1{N}^T{L}^T\mathrm{DLNd\zeta }---(F3-7)$
$P_T^e=l{\int }_0^1{N}^T{L}^{T}D{ϵ}_T\mathrm{d\zeta }---(F3-8)$
如果考虑阻尼和结点力，则上式改写为：
Mee+Ceae+Keae＝Pe                                             (F3-9)
组集之后即可得到总体结构的动力学有限元方程：
M+Ca+Ka＝P                                                    (F3-10)
上式中：
M为结构质量矩阵；
K为结构刚度矩阵；
a为结构节点位移矢量；
C表示阻尼矩阵，通常我们采用比例阻尼(即Rayleigh阻尼)，即：
C＝αM+βK                                                      (F3-11)
其中α、β是不依赖于频率的常数；
P表示结点力列项，对于静坐标系统，它表示外部结点力作用(包括等 效结点力例如重力)，对于动坐标系，它还必须包括由于牵连加速度和柯氏 加速度所产生的牵连力和柯氏力。
下面采用Newmark方法对上述方城在在时间域上离散。Newmark方法实质 上是线性加速度法的一种推广，它采用下列假设：
rt+Δt＝rt+((1-δ)rt+δrt+Δt)                          (F3-12)
$d_{\mathrm{rt}+\mathrm{\Delta t}}=d_{\mathrm{rt}}+{\stackrel{·}{d}}_{\mathrm{rt}}\mathrm{\Delta t}+((\frac{1}{2}-\alpha ){\stackrel{··}{d}}_{\mathrm{rt}}+\alpha {\stackrel{··}{d}}_{\mathrm{rt}+\mathrm{\Delta t}})\Delta t^2---(F3-13)$
其中α和δ是按积分精度和稳定性要求而决定的参数。当α＝1/2且δ＝1/6 时，上面两式相当于线性加速度法，因为这时它们可以从下面时间间隔为Δt内 线性假设的加速度表达式的积分得到：
rt+τ＝rt+(rt+Δt-rt)τ/Δt                              (F3-14)
式中0≤τ≤Δt，Newmark方法原来是从常平均加速度法这样一种无条件稳 定积分方案提出的，那时δ＝1/2且α＝1/4，Δt内的加速度为：
rt+τ＝(rt+rt+Δt)/2                                      (F3-15) Newmark方法中时间t+Δt的位移解答drt+Δt是通过满足t+Δt的运动方程
Mrt+Δt+Cmrt+Δt+Kmdrt+Δt＝Qmt+Δt                         (F3-16)
而得到的。为此，首先从假设中解得：
${\stackrel{··}{d}}_{\mathrm{rt}+\mathrm{\Delta t}}=\frac{1}{\mathrm{\alpha \Delta }{t}^2}(d_{\mathrm{rt}+\mathrm{\Delta t}}-d_{\mathrm{rt}})-\frac{1}{\mathrm{\alpha \Delta t}}{\stackrel{·}{d}}_{\mathrm{rt}}-(\frac{1}{2\alpha }-1){\stackrel{··}{d}}_{\mathrm{rt}}---(F3-17)$
将上式带入运动方程则可得到从drt，rt，rt计算drt+Δt的公式：
$(K_m+\frac{1}{\mathrm{\alpha \Delta }{t}^2}M+\frac{\delta }{\mathrm{\alpha \Delta t}}{C}_m)d_{\mathrm{rt}+\mathrm{\Delta t}}$
$=Q_{\mathrm{mt}+\mathrm{\Delta t}}+m(\frac{1}{\mathrm{\alpha \Delta }{t}^2}{d}_{\mathrm{rt}}+\frac{1}{\mathrm{\alpha \Delta t}}{\stackrel{.}{d}}_{\mathrm{rt}}+(\frac{1}{2\alpha }-1)\stackrel{..}{d_{\mathrm{rt}}})$
$+C_m(\frac{\delta }{\mathrm{\alpha \Delta t}}{d}_{\mathrm{rt}}+(\frac{\delta }{\alpha }-1)\stackrel{·}{d_{\mathrm{rt}}}+(\frac{\delta }{2\alpha }-1)\mathrm{\Delta t}\stackrel{··}{d_{\mathrm{rt}}})---(F3-18)$
需要说明的是，Newmark方法实际上是一种隐式算法，从数学上可以严格 的证明：当δ≥0.5，α≥0.25(0.5+δ)2时，Newmark方法是无条件稳定的，即时间 步长Δt的大小可不影响解的稳定性，但它显然会影响解的精度。
为了方便描述，将上述方程改写成标准有限元形式如下：
$(K^e+\frac{1}{\mathrm{\alpha \Delta }{t}^2}{M}^e+\frac{\delta }{\mathrm{\alpha \Delta t}}{C}^e)a_{t+\mathrm{\Delta t}}^e$
$=Q_{t+\mathrm{\Delta t}}^e+M^e(\frac{1}{\mathrm{\alpha \Delta }{t}^2}{a}_t^e+\frac{1}{\mathrm{\alpha \Delta t}}\stackrel{·}{a_t^e}+(\frac{1}{2\alpha }-1)\stackrel{··}{a_t^e})$
$+C^e(\frac{\delta }{\mathrm{\alpha \Delta }t}{a}_t^e+(\frac{\delta }{\alpha }-1)\stackrel{·}{a_t^e}+(\frac{\delta }{2\alpha }-1)\mathrm{\Delta t}\stackrel{··}{a_t^e})---(F3-19)$
其中：Me、Ce、Ke分别表示单元质量矩阵、单元阻尼矩阵、单元刚度矩 阵，Qe表示单元载荷列阵，ae表示单元结点位移列阵。
注意观察上面的方程，可以仿照静力学方程的形式作进一步改写：
${\stackrel{-}{K}}^e{\stackrel{-}{a}}^e=P^e---(F3-20)$
其中：
${\stackrel{-}{K}}^e=K^e+\frac{1}{\mathrm{\alpha \Delta }{t}^2}{M}^e+\frac{\delta }{\mathrm{\alpha \Delta t}}{C}^e---(F3-21)$
${\stackrel{-}{a}}^e=a_{t+\mathrm{\Delta t}}^e---(F3-22)$
$P^e=Q_{t+\mathrm{\Delta t}}^e+M^e(\frac{1}{\mathrm{\alpha \Delta }{t}^2}{a}_t^e+\frac{1}{\mathrm{\alpha \Delta t}}\stackrel{·}{a_t^e}+(\frac{1}{2\alpha }-1)\stackrel{··}{a_t^e})$
$+C^e(\frac{\delta }{\mathrm{\alpha \Delta t}}{a}_t^e+(\frac{\delta }{\alpha }-1)\stackrel{·}{a_t^e}+(\frac{\delta }{2\alpha }-1)\mathrm{\Delta t}\stackrel{··}{a_t^e})---(F3-23)$
很显然，在按Newmark方法在时间域上离散之后，动力学方程的求解完全 类似于静力学方程。
Newmark方法求解方程(F3-0-5)的方法属于数值计算内容，请使用相关 数值计算方法进行求解，这里不再赘述。
通过求解以上波动方程可以获得各个生产时刻杆柱不同深度处的螺杆泵 工作系统的扭转圈数、速度、扭矩、角度等动态的变化参数，同时杆柱不同 深度处的这些动态参数也对应着摩擦损耗变化，不同深度部分的功率等参数。 也就是说，通过方程求解可以得到如下的扭矩构成和受力构成的各个分量：
M＝M0+M1+M2+M3                                                  (F3-24)
M-光杆扭矩，N.m；M0-螺杆泵初始扭矩，N.m；M1-螺杆泵的工作扭矩， N.m；M2-抽油杆和井液的摩擦扭矩，N.m；M3-抽油杆的惯性扭矩以及扶 正器的摩擦扭矩，N.m。
F＝F1+F2+F3+F4+F5                                               (F3-25)
F1-抽油杆自重，N；F2-抽油杆所受的浮力，N；F3-液体压差作用在转 子上的轴向力，N；F4-液体对抽油杆向上的摩擦力，N；F5-由于过盈引起的 定子对抽油杆的半干摩擦力，N。
需要注意的是，扶正器的摩擦扭矩要依据斜井的井眼轨迹相关计算模型 得到，多相流的计算中也有用到井斜角角度的数据，这些与斜井有关的侧向 摩擦计算与压力梯度计算均用到斜井计算模型其它井类型的斜井部分也全部 用到了该斜井计算模型。斜井计算模型见相关的斜井井眼轨迹计算与侧向力 计算模型，而且石油工程领域这方面的计算模型较多，涉及到斜井的相关计 算，可以任选一计算方法即可，此处不再赘述。
波动方程的求解过程也就是迭代的过程，通过以上叙述可以获得杆柱最 下端的有功功率和排出压力，根据有功功率＝质量流量×排出压力这一基本的 关系式计算出质量流量，这一点和电潜泵井的思想是相同的，所不同的是电 潜泵井油管中只有流体没有杆柱，其泵排出压力使用多相流算出，而螺杆泵 井是在更加复杂的方程求解中逐级求得杆柱受力和多相流压力分布得到有效 功率和排出压力的，也就是说，螺杆泵的有效功率和排出压力是一个统一的 复杂的耦合模型，当然这种耦合是计算机算法上的耦合，其物理计算依据是 相同的。
螺杆泵井得到的是泵工作压力下的地下液体质量流量，折算求出井口有 效排量主要是根据质量连续原理，除以地面条件下的原油密度，就得到地面 体积流量，同时减去气体所占质量分数，一般情况下气体所占质量分数是非 常小的，可以忽略不计。
4.电潜泵井
如图9所示，为电潜泵井产液量计算流程图。根据有功功率＝质量流量× 排出压力这一基本的关系式计算出质量流量。所述工况数据为：瞬时产液量、 累积产液量，输入功率，泵功率、泵效率、有效功率，系统损耗、系统效率、 泵排出压力，电量、电费等，还有丰富的宏观控制图、诊断结论等诸多指标。 所述油井基础数据为嘴前嘴后压力P1，P2、三相电流I1，I2，I3、电压U、功率因素 cosφ等动态参数；油嘴直径d、生产气液比Rs等静态参数根据。
主要根据动态参数：嘴前压力P1(油压)、嘴后压力P2(回压)、三相电 流I1，I2，I3、电压U、功率因数cosφ，静态参数：油嘴直径d、生产气液比Rs； 利用多相流油嘴节流模型，电潜泵、电缆能耗模型，配合举升数学模型加以 修正和拟合，计算出电潜泵井的混合流体流量，再用流量标定系数k计算得到 电潜泵井井口折算体积流量。其理论模型示意函数如下：
Q电潜泵＝kf(d，P1，P2，Rs，I1，I2，I3，U，cosφ)                (F4-0)
Q电潜泵----电潜泵井的产液量，m3；d----油嘴直径，mm；P1----嘴前压力 (油压)，MPa；P2----嘴后压力(回压)，MPa；Ps----生产气液比；I1，I2，I3---- 三相电流，A；U----电压，V；cosφ----功率因数；k----流量标定系数， 小数。
电潜泵井的流量计量基本思想是应用能耗模型获得有功功率，应用多相 流计算出泵的排出压力，应用有功功率＝质量流量×排出压力这一基本的关系 式计算出质量流量，然后根据气油比等参数算出气体所占质量，这个气体质 量比重是非常小的，从而获得液体的质量流量，再根据区块原油的密度可以 获得该质量流量原油对应的地面液体体积流量，也就是产液量。电潜泵井得 到的是泵工作压力下的地下液体质量流量，这一点和螺杆泵相同。
电潜泵井的流量计量获得有功功率的能耗模型如以下模型所示，获得泵 的排出压力的多相流计算模型与自喷井中的模型一致，此处不再赘述。
有效举升高度(m)：H有效＝H中+(P油-P流)×102                      (F4-1)
有效功率(kw)：P有效＝Q×H有效×γ混×100/8812800                (F4-2)
有效扬程(m)：H有效扬程＝H有效+HL                                (F4-3)
油管压头损失(MPa)：HL＝0.111×10-10×λ×H泵深×q2/d5，         (F4-4)
其中，q为平均流体体积流量，q＝Q/γ混，注意应该γ混＝0.5---1之间，λ可 以取常数0.06；
油管摩阻损失功率(kw)：P油管＝P有效×HL/H有效                    (F4-5)
泵吸入口压力(MPa)：P泵入口＝P流-γ混(H中-H泵)×0.00981          (F4-6)
电缆压降(V)：
电机输入电压(V)：U电机入＝U-ΔU电缆                             (F4-8)
电机输入功率(kw)：电机输入功率＝电缆输入功率-电缆消耗功率
电缆输入功率(kw)：P输入＝1.732×U×(I1+I2+I3)cosφ/3            (F4-9)
电缆消耗功率(kw)：
L为电缆长度(m)，
其中，R电缆＝ρ/S                                               (F4-11)
r＝k*p/s (温度＝(t口+t底)/2下的k值)
其中，ρ为铜导线电阻率(0.30295)，S为导体截面积，
地面设备损失功率(kw)：P地面＝P输入×0.025                       (F4-12)
日耗电量(度)：W日耗＝P输入×1.025×24                           (F4-13)
方液日耗电量(度)：W方液日耗＝Ｗ日耗/Q                           (F4-14)
系统效率(％)：η＝P有效×100％/P输入                            (F4-15) 混合液密度由下式求得：

其中， fw为含水，fw＝φ/100。
通过本发明的系统和方法，可完成的功能为：
1.油井工况检测功能
采集自喷井、抽油机井、电潜泵井、螺杆泵井电压、电流、功率、载荷、 冲次、冲程、井口压力、油温、生产时率、曲柄销子、转速、扭矩及巡井时 间等生产参数，并实现人工/自动远程控制。
2.注水井工况检测功能
对注水井采集井口压力、温度、注水量、巡井时间等生产参数，并实现 人工/自动远程控制。
3.故障报警功能
停电、停机、回压异常、缺相及电流异常、抽油机抽空、防盗红外监控， 曲柄销松动脱落。
4.控制保护功能
对油井间抽控制；缺相、三相电流不平衡、曲柄销脱落自动停机；远程 启停控制和开关井场照明灯；实现对油水井产/注液量的自动计量控制。
5.数据通信功能
单井采用自主开发的YDSW远程数据采集控制器(RTU)与上位机进行数据 通信；中控室采用无线宽带/光缆网络通讯方式实现联网。
6.数据管理功能
图形方式实时显示压力、温度、载荷、扭矩、电流、电压、功率等各种 生产参数及泵、机等生产设备运行状态；实现生产参数超限报警及设备故障 报警，预测故障位置和故障原因并进行相应提示。
7.摄像监视功能
在油水井、接转站或联合站外安装若干云台和全景低照度摄像机，对站 内全貌及油水井进行监视。
8.生产管理及遥控指挥
自动记录巡井时间；与油田局域网数据共享；可以通过现有局域网络， 在网上远端监控油水井生产现场并进行指挥。
9.油井产量计量功能
在油井正常运行期间，能够实现远程测试数据的自动录取，根据采集的 压力、转速、地面示功图等数据，应用油井量油技术，计算油井产液量，实 现在无人值守情况下能及时掌握油井的动态变化。
10.油井生产系统分析与优化决策功能
根据检测数据，进行生产井参数优化设计、在线诊断、抽油井系统效率 分析等。
11.网络查询功能
通过IE浏览器和专用的视频播放软件，在油田信息网上可随时浏览各油 水井的各个监控画面及实时生产数据，及液量计算结果，查询有关生产报表 及分析结果。
通过本发明，可在油田信息网上全面了解工区的生产状况，实现油水井 的远程实时监控、联合站生产参数和设备运行状态的自动检测与重要生产岗 位的闭环控制和图像监视，提高油田生产时率，简化地面流程、节能降耗， 保证油田安全生产，提高经济、社会效益。
上述实施例仅用于说明本发明，而非用于限定本发明。