MATLAB实现有杆抽油系统的数学建模及诊断_gibbs模型将悬点示功图转化为泵功图-CSDN博客

本文链接：https://blog.csdn.net/sheziqiong/article/details/136428557

本文详细探讨了光杆悬点的运动规律，通过数学模型推导得到位移、速度和加速度函数，并应用Gibbs模型将悬点示功图转化为泵功图。此外，文中还介绍了如何利用泵功图计算一级杆油井的产量，涉及有效冲程法和气体检测模型。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

一、问题重述 3
二、问题假设 4
三、符号说明 5
四、问题一：光杆悬点运动规律 5
4.1问题分析 5
4.2模型建立及推导 6
4.2.1悬点运动参数方程 6
4.2.2悬点运动速度 7
4.2.3悬点运动加速度 8
4.2.4模型初始条件 8
4.2.5阈值范围及初始条件 8
4.2.6模型验证及误差分析 10
4.3 结论 12
五、问题二：泵功图计算 12
5.1Gibbs 模型的基本思想 12
5.2Gibbs 波动方程的推导 12
5.3Gibbs 波动方程的求解 14
5.3.1一级杆 Gibbs 波动方程的求解 14
5.3.2多级杆 Gibbs 波动方程的求解 16
5.4模型详细实现过程 17
5.4.1原始数据的处理 17
5.4.2边界条件 18
5.4.3初始条件 18
5.4.4核心计算步骤 18
5.4.5程序设计流程图 20
5.5模型的验证 21
5.6一级杆油井的模型验证 21
5.7三级杆油井的模型验证 22
5.8 小结 22
六、问题三：泵功图的应用 23
6.1背景分析与假设 23
6.2问题一有效冲程法模型 23
6.2.1模型的建立与推导 23
6.2.2模型的验证 27
6.2.3产量计算 31
6.3问题一有效功率求解法模型 32
6.3.1模型的建立与推导 32
6.3.2模型的验证 36
6.3.3产量计算 36
6.4问题二模型 36
6.4.1模型的推导 36
6.4.2模型的建立 37
6.4.5 模型使用 37
七、问题四：Gibbs 模型深入研究 40
7.1问题二 40
7.1.1基于摩擦功的阻尼系数模型建立 40
7.1.2模型求解 43
八、模型评价 44
九、参考文献 45
一、问题重述

目前，开采原油广泛使用的是有杆抽油系统。电机旋转运动转化为抽油杆上下往返周期运动，带动设置在杆下端的泵的两个阀的相继开闭，从而将地下上千米深处蕴藏的原油抽到地面上来。
在一个冲程期间，仪器以一系列固定的时间间隔测得悬点处的一系列位移数据和荷载数据，据此建立悬点的示功图称为悬点示功图。因为受到诸多因素的影响，在同一时刻，悬点处的荷载与柱塞的受力以及悬点、柱塞处的相对位移都不相同，因此悬点示功图与泵功图是不同的。工程上一般根据示功图形状与理论示功图进行对比来判断抽油机工作状态。通过悬点示功图可以初步诊断该井的工作状况，如产量、气体影响、等。要精确诊断油井的工作状况，需要采用泵功图。通过数学建模，把悬点示功图转化为杆上任意点的示功图，并最终确定泵功图。

问题一：光杆悬点运动规律
电机旋转运动通过四连杆机构转变为抽油杆的垂直运动。悬点下只挂光杆。根据已给的四连杆各段尺寸，利用附件参数，求出悬点的一个冲程的运动规律：位移函数、速度函数、加速度函数。并与有荷载的附件悬点位移数据进行比较。问题二：泵功图计算
使用 Gibbs 模型，给出由悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；根据附件 1、3 油井参数和悬点示功数据，分别计算这两口油井的泵功图数据，并绘制出两油井的悬点示功图和泵功图。
问题三：泵功图的应用
1）建立 2 个不同的由泵功图估计油井产量的模型，至少一个要利用“有效冲程”；并利用附件数据估算两口油井一天（24 小时）的产液量。
2）当泵中有气体，导致泵没充满。建立模型或算法，自动判别某泵功图数据是否属于泵内有气体的情况。并对附件对应的泵功图进行计算机诊断是否属于泵内充气这种情况。

%-----------------------------------------------------------
%题号:  问题二—泵功图计算之一级杆某油井的泵功图计算
%参考数据： 附件1  7#井
%作者:  102686020  陈源源，游检卫，王旭阳
%时间： 2012.09.21-2012.09.24 
%Email: jvyou@seu.edu.cn
%----------------------------------------------------------

%%
clear all;
clc;

%%
%-----加载外部数据（边界条件）-------
data1=load('1.txt');
u_1=data1(1,:);   %悬点处位移时间函数；
F_1=data1(2,:)*1e3;   %悬点处荷载时间函数；
% figure(1)
% plot(u_1,F_1,'r*-');grid on;
% xlabel('位移 (m)');ylabel('荷载 (n)');title('一级杆的悬点示功图');
% legend('悬点示功图');
rod_Num=1;

%%
%---------简单参数的初始化------------
data_length=length(u_1);  %K=K1=data_length

E=2.1*1e11;
density=8456;
a_velocity=sqrt(E/density);

radius=(22*(1e-3))/2;
A=pi*radius^2;

n_ba_u=20;
n_x=2;      n_t=data_length;

u=zeros(n_x,n_t);F=zeros(n_x,n_t);
u(1,:)=u_1;  F(1,:)=F_1;
L=792.5;     T=(1/7.6)*60;
delta_t=T/(n_t-1);
omega=2*pi/T;

%%
% %----[1]【a】基于Gibbs模型的阻尼系数c计算---
% S=3.2;
% N=7.6;
% V=S*N/30;
% % gama=exp(0.18-0.39094*V*30);
% gama=0.19746-0.06974*V+0.00075*V^2;
% c=(pi*a_velocity*gama)/(2*L);

% % %----[1]【b】基于摩擦功的阻尼系数c计算---
Mur_fluid=30*1e-3;  %液体动力粘度
L_crank=950*1e-3;   %驱柄杆长度
Dt=70*1e-3;         %油管内径
Dr=22*1e-3;         %抽油杆直径
m_Dia=Dt/Dr;        %直径系数
Hp=793;         %泵深
B1=(m_Dia^2-1)/(2*log(m_Dia))-1;
B2=m_Dia^4-1-(m_Dia^2-1)^2/log(m_Dia);
% % % 方法一
% c=(2*pi*Mur_fluid/(9.8*density*A))*(1/log(m_Dia)+(4*(B1+1)/B2)*...
%     (B1+(2*tan(omega*L_crank/a_velocity))/((omega*L_crank/a_velocity)/...
%     cos((omega*L_crank/a_velocity))+sin((omega*L_crank/a_velocity)))));
% % 方法二
c=(2*pi*Mur_fluid/(density*A))*(1/log(m_Dia)+(4*(B1+1)/B2)*...
    (B1+(2*tan(omega*Hp/a_velocity))/((omega*Hp/a_velocity)/...
    sin((omega*Hp/a_velocity))+cos((omega*Hp/a_velocity)))));
%%
alpha_n=zeros(1,n_ba_u);
belta_n=zeros(1,n_ba_u);
%----alpha_n,belta_n的计算----
for iN=1:n_ba_u+1
    alpha_n(1,iN)=(iN*omega/(a_velocity*sqrt(2)))*sqrt(1+sqrt(1+(c/(iN*omega))^2));
    belta_n(1,iN)=(iN*omega/(a_velocity*sqrt(2)))*sqrt(-1+sqrt(1+(c/(iN*omega))^2)); 
end

%%
%--------[2]系数的计算--------
%--------sigma_n,tao_n,v_n,delta_n的计算--------
sigma_n=zeros(1,n_ba_u+1);
tao_n=zeros(1,n_ba_u);
v_n=zeros(1,n_ba_u+1);
delta_n=zeros(1,n_ba_u);
for iN=1:n_ba_u+1  
    %%
    sigma_n_iNsum=0; v_n_iNsum=0;
    %-----sigma_n,v_n的计算----
    for jP=1:data_length
        sigma_n_iNsum=sigma_n_iNsum+F_1(jP)*cos(2*(iN-1)*pi*jP/data_length);
        v_n_iNsum=v_n_iNsum+u_1(jP)*cos(2*(iN-1)*pi*jP/data_length);
    end
%     sigma_n(1,1)=sigma_n(1,1)/data_length;
    sigma_n(1,iN)=(2/data_length.^2)*sigma_n_iNsum;
    v_n(1,iN)=(2/data_length)*v_n_iNsum;
    
    %%
    if iN<=n_ba_u
       tao_n_iNsum=0; delta_n_iNsum=0;
       %-----tao_n,delta_n的计算----   
       for jP=1:data_length
          tao_n_iNsum=tao_n_iNsum+F_1(jP)*sin(2*(iN)*pi*jP/data_length);
          delta_n_iNsum=delta_n_iNsum+u_1(jP)*sin(2*(iN)*pi*jP/data_length);
       end
       tao_n(1,iN)=(2/data_length)*tao_n_iNsum;
       delta_n(1,iN)=(2/data_length)*delta_n_iNsum;
    end        
end 

%%
%-----[3]Kai_n,Mur_n的计算----
Kai_n=zeros(1,n_ba_u);
Mur_n=zeros(1,n_ba_u);
for iN=1:n_ba_u
    Kai_n(1,iN)=(sigma_n(1,iN+1)*alpha_n(1,iN)+tao_n(1,iN)*belta_n(1,iN))/(E*A*(alpha_n(1,iN)^2+belta_n(1,iN)^2));
    Mur_n(1,iN)=(sigma_n(1,iN+1)*belta_n(1,iN)-tao_n(1,iN)*alpha_n(1,iN))/(E*A*(alpha_n(1,iN)^2+belta_n(1,iN)^2));    
end

%%
%---------[4]核心计算---------------
%%
%-----[4.1]----泵处--On,Pn,On_1,Pn_1的计算------
On=zeros(1,n_ba_u);
Pn=zeros(1,n_ba_u);
On_1=zeros(1,n_ba_u);
Pn_1=zeros(1,n_ba_u);
for iN=1:n_ba_u
    %%
    alpha_n_iN=alpha_n(1,iN);    belta_n_iN=belta_n(1,iN);
	sigma_n_iN=sigma_n(1,iN+1);  tao_n_iN=tao_n(1,iN);      v_n_iN=v_n(1,iN+1);     delta_n_iN=delta_n(1,iN);
	Kai_n_iN=Kai_n(1,iN);        Mur_n_iN=Mur_n(1,iN);
    %%
    On(1,iN)=(Kai_n_iN*cosh(belta_n_iN*L)+delta_n_iN*sinh(belta_n_iN*L))*sin(alpha_n_iN*L)...
                +(Mur_n_iN*sinh(belta_n_iN*L)+v_n_iN*cosh(belta_n_iN*L))*cos(alpha_n_iN*L);
    Pn(1,iN)=(Kai_n_iN*sinh(belta_n_iN*L)+delta_n_iN*cosh(belta_n_iN*L))*cos(alpha_n_iN*L)...
                -(Mur_n_iN*cosh(belta_n_iN*L)+v_n_iN*sinh(belta_n_iN*L))*sin(alpha_n_iN*L);  
    On_1(1,iN)=((delta_n_iN*belta_n_iN-v_n_iN*alpha_n_iN)*cosh(belta_n_iN*L)+tao_n_iN*sinh(belta_n_iN*L)/(E*A))*sin(alpha_n_iN*L)...
            +((v_n_iN*belta_n_iN+delta_n_iN*alpha_n_iN)*sinh(belta_n_iN*L)+sigma_n_iN*cosh(belta_n_iN*L)/(E*A))*cos(alpha_n_iN*L);
    Pn_1(1,iN)=((delta_n_iN*belta_n_iN-v_n_iN*alpha_n_iN)*sinh(belta_n_iN*L)+tao_n_iN*cosh(belta_n_iN*L)/(E*A))*cos(alpha_n_iN*L)...
            -((v_n_iN*belta_n_iN+delta_n_iN*alpha_n_iN)*cosh(belta_n_iN*L)+sigma_n_iN*sinh(belta_n_iN*L)/(E*A))*sin(alpha_n_iN*L);   
end

%%
%-----[4.2]----泵处--u(L,t),F(L,t)的计算------
for nt=1:n_t
    %%
    sum_u=0;    sum_F=0;
    for iN=1:n_ba_u
        sum_u=sum_u+On(1,iN)*cos(iN*omega*(nt-1)*delta_t)+Pn(1,iN)*sin(iN*omega*(nt-1)*delta_t);
        sum_F=sum_F+On_1(1,iN)*cos(iN*omega*(nt-1)*delta_t)+Pn_1(1,iN)*sin(iN*omega*(nt-1)*delta_t);
    end        
     %%
    sum_u_t(nt)=sum_u;
    sum_F_t(nt)=sum_F;
    EA=E*A;
    %--位移时间函数和荷载时间函数的计算----
    %位移时间函数的计算u(x,t)
    u(2,nt)=sigma_n(1,1)*L/(2*E*A)+v_n(1,1)/2+sum_u;
    %荷载时间函数的计算F(x,t)
    F(2,nt)=sigma_n(1,1)/2+E*A*sum_F;
end

%%
%---------示功图的显示------------
figure(2)
plot(u(1,:),F(1,:),'k-');
hold on;

plot(u(2,:),F(2,:),'r-.');
grid on;
xlabel('位移 (m)');ylabel('荷载 (n)');title('一级杆的悬点示功图');
legend('悬点示功图','泵功图');
% hold off

在这里插入图片描述