MATLAB实现有杆抽油系统的数学建模及诊断1

目 录

1问题重述 6
2问题分析 7
2.1问题一：光杆悬点运动规律 7
2.2问题二：泵功图计算 7
2.3问题三：泵功图的应用 8
2.4问题四：深入研究的问题 8
3符号定义与说明 9
4问题假设 11
5问题一：光杆悬点运动规律 12
5.1简化模型 12
5.1.1简化为简谐运动模型 12
5.1.2简化为曲柄滑块机构模型 13
5.2解析模型 15
5.3计算结果比较及分析 17
6问题二：泵供图计算 19
6.1计算模型及求解方法 19
6.2计算结果与分析 22
7问题三：泵功图的应用 24
7.1油井日产液量计算 24
7.1.1模型一：有效冲程法 24
7.1.2模型二：面积法 30
7.1.3计算结果及分析 32
7.2泵故障诊断 32
7.2.1泵功图网格化 32
7.2.2灰度矩阵的形成 33
7.2.3灰度矩阵特征值的提取 33
7.2.4建立灰关联度 34
7.2.5计算结果及分析 35
8问题四：深入研究的问题 36
8.1抽油系统模型改进 36
8.1.1非齐次波动方程模型建立及求解 36
8.1.2计算结果与分析 38
9模型评价与改进 41
9.1问题一模型评价 41
9.2问题二模型评价 41
9.3问题三模型评价 41
9.4问题四模型评价 41
10参考文献 42
1 问题重述

问题一：光杆悬点运动规律
电机旋转运动通过四连杆机构转变为抽油杆的垂直运动。假设驴头外轮廓线为部分圆弧、电机匀速运动，悬点 E 下只挂光杆（光杆下不接其它杆，不抽油，通常用来调试设备）。请按附录 4 给出四连杆各段尺寸，利用附件 1 的参数，求出悬点 E 的一个冲程的运动规律：位移函数、速度函数、加速度函数。并与有荷载的附件 1 的悬点位移数据进行比较。
问题二：泵功图计算
1966 年，Gibbs 给出了悬点示功图转化为地下示功图的模型，由于受计算机速度的限制，直到近些年才得以被重新重视。请使用 Gibbs 模型，给出由悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；附件
1 是只有一级杆的某油井参数和悬点示功数据，附件 2 是有三级杆的另一油井参数和悬点示功数据，利用它们分别计算出这两口油井的泵功图数据；并分别绘制出两油井的悬点示功图和泵功图（每口井绘一张图，同一井的悬点示功图与泵功图绘在同一张图上， 请标明坐标数据）。
问题三：泵功图的应用
1）建立 2 个不同的由泵功图估计油井产量的模型，其中至少一个要利用“有效冲程”；并利用附件 1 和附件 2 的数据分别估算两口油井一天（24 小时）的产液量。（单位：吨， 这里所指的液体是指从井里抽出来的混合液体）
2）如图 5（C）形式的泵功图表示泵内有气体，导致泵没充满。请建立模型或算法， 以由计算机自动判别某泵功图数据是否属于泵内有气体的情况。并对附件 1、附件 2 对应的泵功图进行计算机诊断是否属于泵内充气这种情况。
问题四：深入研究的问题
1）请对 Gibbs 模型进行原理分析，发现它的不足。在合理的假设下，重新建立抽油系统模型或对现有模型进行改进；并给出由悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；利用附件 1、附件 2 的数据重新进行计算；对计算结果与问题二的计算结果进行比较，分析你的模型的优缺点。
2u 2 2u u

2）Gibbs 模型在数学上可简化为 “波动方程”： t 2  a

x2  c t

其中 a 为已知

常数，c 称为阻尼系数，鉴于大多数的阻尼系数公式是作了诸多假设后推出的，并不能完整地反应实际情况。如果能从方程本身和某些数据出发用数学方法估计参数 c，贡献是很大的。对此，请你进行研究，详细给出计算 c 的理论推导过程并尽可能求出 c。如果需要题目之外的数据，请用字母表示之并给出计算 c 的推导过程。

function Code21
clc;
%xx=8456*(pi*22*22/4/1000/1000)*792.5
A=load('test1.txt');
U=A(:,1);
D=A(:,2)*1000; %原始载荷
k=length(D);
n=10;
plot(U,D,'r-')
hold on

%计算声速，粘滞阻力系数
EE=2.1e11;
ru=8456;
AA=sqrt(EE/ru);
DL=70/1000.0;
DR=22/1000.0;
LL=792.5;        %一级杆杆长(m)
mm=30/1000.0;    %地面原油粘度(Pa*s)
cc=7.6;          %冲数
ww=2*pi*cc/60.0; %转速

rul=864*(1-0.98)+1000*0.98;
D=D-(ru-rul)*(pi*DR*DR/4)*LL*9.8; %去重力载荷
%计算sig0,vel0
sig0=0.0;
vel0=0.0;
for i=1:k
    sig0=sig0+2/k*D(i);
    vel0=vel0+2/k*U(i);
end
%sig(i),tau(i),vel(i),det(i)
for i=1:n
    sig(i)=0.0;
    tau(i)=0.0;
    vel(i)=0.0;
    det(i)=0.0;
    for p=1:k
       sig(i)=sig(i)+2/k*D(p)*cos(2*pi*i/k*p); 
       tau(i)=tau(i)+2/k*D(p)*sin(2*pi*i/k*p);
       vel(i)=vel(i)+2/k*U(p)*cos(2*pi*i/k*p); 
       det(i)=det(i)+2/k*U(p)*sin(2*pi*i/k*p);
    end
end

%test
for p=1:k
    DD(p)=sig0/2;
    UU(p)=vel0/2;
    for i=1:n
       DD(p)=DD(p)+sig(i)*cos(p*i/k*2*pi)+tau(i)*sin(p*i/k*2*pi);
       UU(p)=UU(p)+vel(i)*cos(p*i/k*2*pi)+det(i)*sin(p*i/k*2*pi);
    end
end
%plot(UU,DD,'+')

%粘滞阻力系数
DR=22/1000.0;
LL=792.5;        %一级杆杆长(m)
pm=DL/DR;
B1=(pm^2-1)/2/log(pm)-1;
B2=pm^4-1-(pm^2-1)^2/log(pm);
tp=ww*LL/AA*sin(ww*LL/AA)+cos(ww*LL/AA);
pc=8*mm/ru/DR/DR*(1/log(pm)+2/B2*(B1+1)*(B1+2/(tp)));
%计算ahl(i),bet(i),kin(i),muu(i)
EA=EE*(pi*DR*DR/4);
for i=1:n
    aph(i)=i*ww/AA/sqrt(2)*sqrt(1+sqrt(1+(pc/i/ww)^2));
    bet(i)=i*ww/AA/sqrt(2)*sqrt(-1+sqrt(1+(pc/i/ww)^2));
    kin(i)=(sig(i)*aph(i)+tau(i)*bet(i))/EA/(aph(i)^2+bet(i)^2);
    muu(i)=(sig(i)*bet(i)-tau(i)*aph(i))/EA/(aph(i)^2+bet(i)^2);
end
%计算OO,PP,DO,DP
LL=793; %泵深(m)
for i=1:n
    OO(i)=(kin(i)*cosh(bet(i)*LL)+det(i)*sinh(bet(i)*LL))*sin(aph(i)*LL)...
       +(muu(i)*sinh(bet(i)*LL)+vel(i)*cosh(bet(i)*LL))*cos(aph(i)*LL); 
    PP(i)=(kin(i)*sinh(bet(i)*LL)+det(i)*cosh(bet(i)*LL))*cos(aph(i)*LL)...
       -(muu(i)*cosh(bet(i)*LL)+vel(i)*sinh(bet(i)*LL))*sin(aph(i)*LL);
    tp1=tau(i)/EA; %bet(i)*kin(i)-aph(i)*muu(i); %tau(i)/EA;
    tp2=det(i)*bet(i)-vel(i)*aph(i);
    tp3=sig(i)/EA; %aph(i)*kin(i)+bet(i)*muu(i); %sig(i)/EA;
    tp4=vel(i)*bet(i)+det(i)*aph(i);
    DO(i)=(tp1*sinh(bet(i)*LL)+tp2*cosh(bet(i)*LL))*sin(aph(i)*LL)...
       +(tp3*cosh(bet(i)*LL)+tp4*sinh(bet(i)*LL))*cos(aph(i)*LL);
    DP(i)=(tp1*cosh(bet(i)*LL)+tp2*sinh(bet(i)*LL))*cos(aph(i)*LL)...
       -(tp3*sinh(bet(i)*LL)+tp4*cosh(bet(i)*LL))*sin(aph(i)*LL);
end
%计算位移
tp0=sig0/2/EA;
for p=1:k
   Xx(p)=p/k*360;
   Ux(p)=tp0*LL+vel0/2;
   Fx(p)=sig0/2;
   for i=1:n
      Ux(p)=Ux(p)+OO(i)*cos(p*i/k*2*pi)+PP(i)*sin(p*i/k*2*pi);
      Fx(p)=Fx(p)+EA*(DO(i)*cos(p*i/k*2*pi)+DP(i)*sin(p*i/k*2*pi));
   end
end
%plot(Xx,Ux,'+')

%减油压，套压差
%Fp=(0.3-0.2)*pi/4*70^2;
%Fx=Fx-Fp;

plot(Ux,Fx,'g*')
hold on
fid = fopen('case21.txt', 'wt');
for i=1:k
    fprintf(fid, '%f\t%f\n', Ux(i),Fx(i));
end 
fclose(fid);