目录
题 目 有杆抽油系统的数学建模及诊断 1
摘 要 1
— 问题的重述与分析 3
1.1 问题重述 3
1.2 需要解决的问题 3
二 问题假设 4
三 符号说明 5
四 模型建立与求解 5
4.1 问题一：光杆悬点的运动规律 5
4.1.1 模型建立 5
4.1.2 简化为曲柄滑块机构求解悬点 E 运动规律 4
则 5
4.1.3 由解析式求得运动规律的精确解 6
s2  b2  l 2 ds 8
4.1.4 两种方法结果对比 9
12 时悬点到达下死点，下死点位移为 9
4.1.5 对问题一的再思考 10
12 时悬 11
4.2 问题二泵功图计算 11
4.2.1 Gibbs 模型介绍 11
2Efr 12
2 K 2 2n 12
2 K 2 2n 12
2 K 2 2n 12
2 K 2 2n 12
x  (     ) sinh  x) cos x 13
 14
4.2.2 Gibbs 模型在计算泵功图中的应用 14
4.2.2.1 Gibbs 模型理论应用 14
4.2.2.2 利用 Gibbs 模型求解步骤 14
4.2.3 模型求解结果及分析 17
4.2.3.1 对附件 1 中 1#抽油井求解泵功图 17
4.2.3.2 对附件 2 中 7#号抽油井求解泵功 19
2.5 21
4.3 问题三 估计产油量 21
4.3.1 模型 1：“面积法”估计油井产量 21
4.3.1.1 模型分析 21
4.3.2 模型 2：有效冲程估算油井产量 24
4.3.2.1 泵功图模型再分析 24
4.3.2.2 曲率模型的建立 25
i L L L L L L 26
4.3.2.3 柱塞有效冲程的确定方法 26
4.3.2.4 抽油机井产液量计算 27
= 1-fw  28
4.3.3 模型仿真与分析 28
3.5 28
4.4 问题四 深入研究的问题 30
4 Gibbs模型求解泵功图(c=2.5) 31
4.4.1 Gibbs 模型的改进 31
4.4.2 阻尼系数计算方法 35
2L 35
(PRhp  Hhp)(T / S )2 35
4.4.2.2 基于示功图的阻尼系数计算方法 35
u(x, t) dx 35
0 1 AP 39
4.4.2.3 实验仿真 40
五 模型的评价与推广 41
— 问题的重述与分析

1.1问题重述
有杆泵往复抽汲方式是应用最早也最为广泛的一种人工举升采油法，其具有设备装置简单、操作方便、综合成本低的特点。早在石油工业问世时，就开始采用这一方法进行采油。目前，在各种人工举升采油方法中，有杆抽油仍居于首要地位。然而，抽油井井下工作环境复杂，易于发生故障，且井下难于长期放置监测装置，使得井下故障难于及时发现， 因此开发一套抽油井诊断系统，对油田降低能耗，提高生产效率具有重要意义。
1966 年美国壳牌石油公司的S.G.Gibbs 和A.B.Neely 提出了有杆抽油系统的计算机诊断技术，并于 1967 年 9 月 26 日获得美国专利权[1][2]。其基本原理是把抽油杆柱看作井下动态信号的传导线，而抽油泵则为发送器，泵的工作状况以应力波的形式沿抽油杆柱传递到地面，被作为接受器的动力仪所接收。根据这个原理，建立了带阻尼的波动方程作为描述应力波在抽油杆柱中传递过程的基本微分方程，以光杆动载荷及位移作为边界条件求解此方程，可以得到抽油杆柱任意截面处的位移和载荷，从而绘出所需的抽油杆柱任意截面及泵的示功图。根据这些示功图，可以对抽油系统的工作状况作出分析和判断。
1.2需要解决的问题
1、建立杆式抽油机数学模型，求出悬点 E 的一个冲程的运动规律：位移函数、速度函数、加速度函数。
2、利用 Gibbs 模型，结合给出的附件数据和油井参数，通过对数据的预处理，明确边界， 将已有的悬点示功图转化为泵示功图，并画出图形。
3、建立不同的由泵功图估计油井产量的模型，并利用模型估算附件中两口油井一天的产液量。
4、分析 Gibbs 模型的优缺点，尝试体会出新的模型，根据已有数据计算泵功图，并将其与
Gibbs 模型进行对比；目前由于阻尼系数的求解多在假设基础之上，不能完整的反应实际情况，尝试从 Gibbs 波动方程本身或者是现有数据给出阻尼系数的计算方法。

clc;
clear;

% e=2.1e+11;
% load uu1;
% load dd1;
% cc=7.6; %冲次
% x=[792.5];  %抽油杆杆长

load uu2;
load dd2;
cc=4; %冲次
x=[523.61 664.32 618.35];

uu=-uu;
len=length(uu);
pt=find(uu==min(uu));
point=pt(length(pt));
uu=[uu(point:len) uu(1:point-1)];
dd=[dd(point:len) dd(1:point-1)];

ss=0;
l1=0;
uu0=uu(1);
dd0=dd(1);
m=2;
n=length(uu);

while(1)
    deltam=uu(m)-uu0;
    deltan=uu(n)-uu0;
    if deltam<=deltan
        l2=dd(m)-dd0+(dd0-dd(n))/deltan*deltam;
        uu0=uu(m);
        dd0=dd(m);
        delta=deltam;
        m=m+1;
    end
    if deltam>deltan
        l2=dd0-dd(n)+(dd(m)-dd0)/deltam*deltan;
        uu0=uu(n);
        dd0=dd(n);
        delta=deltan;
        n=n-1;
    end
    if(m>n)
        break;
    end
    ss=ss+(l1+l2)*delta/2;
    l1=l2;
end

product=1.44*cc*ss/(sum(x)*9.8);