基于功能原理的有杆抽油系统的数学建模及诊断

目 录
1问题重述 4
1.1研究背景 4
1.2问题一：光杆悬点运动规律 4
1.3问题二：泵功图计算 4
1.4问题三：泵功图的应用 4
1.5问题四：深入研究的问题 5
2问题分析 5
3模型的假设及符号说明 6
3.1模型假设 6
3.2符号说明 6
4光杆悬点运动规律 7
4.1简化为简谐运动时悬点运动规律 7
4.2简化为曲柄滑块机构时悬点运动规律 7
4.3悬点实际运动规律 9
5泵功图的计算与绘制 13
5.1基于 Gibbs 模型的泵功图的计算 13
5.2示功图与泵功图的绘制 18
6泵功图的应用 22
6.1泵功图确立固定阀和游动阀开闭临界点 22
6.1.1泵功图的预处理 22
6.1.2固定阀和游动阀开闭临界点的确立 22
6.2利用泵功图估计油井产量 24
6.2.1功转化分析估算油井产量 25
6.2.2利用有效冲程估算油井产量 26
6.3利用泵功图诊断泵内是否含有气体 27
6.3.1气体对泵功图的影响分析 27
6.3.2面积判别法 27
7问题的进一步研究 28
7.1Gibbs 模型的改进 28
7.1.1基于重力分布载荷的改进 Gibbs 模型 29
7.1.2有限差分法求解 29
7.2基于能量守恒的阻尼系数模型与求解 33
8模型的评价 36
参考文献 37
1 问题重述

1.1研究背景

在原油开采中，有杆抽油系统被广泛使用。电机旋转运动转化为抽油杆上下往返周期运动，带动设置在杆下端的泵的两个阀的相继开闭，从而将地下上千米深处蕴藏的原油抽到地面上来。
钢制抽油杆有一级杆和多级杆，描述抽油杆中任意一水平截面处基本信息的通用方法是示功图，其函数关系表现为位移-荷载关于时间的参数方程。抽油杆上端点称为悬点，在一个冲程期间，仪器以一系列固定的时间间隔测得悬点处的一系列位移数据和荷载数据，据此建立悬点的示功图称为悬点示功图。
“泵”是由柱塞、游动阀、固定阀、部分油管等几个部件构成的抽象概念， 泵中柱塞处的示功图称为泵功图。由于受到诸多因素的影响，在同一时刻t ，悬点处的荷载及相对位移与柱塞的都是不相同的，因此悬点示功图与泵功图是不同的。示功图包含了很多信息，可通过示功图特征分析判断油井的工作状态。
通过悬点示功图可以初步诊断该井的工作状况，如产量、气体影响、阀门漏液、沙堵等，而泵功图能更精确诊断油井的工作状况。然而，泵在地下深处，使用仪器测试其示功数据实现困难大、成本高。因此，建立模型，把悬点示功图转化为杆上任意点的示功图并最终确定泵功图，以准确诊断该井的工作状况，是一个很有价值的实际问题。

clear;

%% 参数
p_water = 0.98;                     % 含水率
chongci = 7.6;                      % 冲次
x = [792.5/3 792.5/3 792.5/3 ];     % 各节长
R = [22 22 22]/1000/2;              % 各节截面半径
gamma = 0.2304;                     % 阻尼系数(无量纲)

% 一个冲程中, 悬点载荷、悬点位移的值(等时间间隔采样)
L = [39.7,39.665,39.665,38.131,31.822,30.184,28.372,27.047,26.385,23.806,19.17,16.765,14.988,13.245,11.885,14.221,17.114,19.693,20.808,22.446,23.422,21.087,18.996,16.451,15.545,13.698,13.001,14.465,16.277,18.159,17.183,17.323,17.288,16.521,15.371,14.709,15.057,14.43,14.256,14.36,14.534,14.709,14.883,14.778,14.813,14.709,14.953,15.266,15.406,15.301,15.51,16.242,16.73,17.253,17.671,18.612,19.031,19.344,20.146,20.808,21.784,21.959,22.69,23.213,23.492,23.353,23.388,23.806,24.12,24.538,24.886,25.339,28.546,31.16,33.251,35.412,37.957,40.78,46.148,50.4,53.99,55.942,56.883,57.092,57.545,58.103,58.138,58.242,58.347,58.068,58.068,58.242,58.661,58.905,59.462,60.02,60.891,61.205,61.344,60.891,61.275,61.066,60.752,59.811,59.497,58.87,58.033,57.615,57.232,56.465,55.524,55.036,54.408,53.955,53.502,53.537,53.049,53.223,52.7,52.247,52.178,52.247,51.481,50.783,50.33,49.807,49.041,48.274,47.298,46.74,45.834,45.59,45.207,44.963,44.684,44.405,44.196,43.917,43.603,43.115,42.732,41.965,40.85,40.257
] .* 1000;
U = [2.518,2.512,2.504,2.495,2.484,2.472,2.459,2.444,2.427,2.41,2.391,2.372,2.351,2.329,2.305,2.281,2.256,2.229,2.201,2.172,2.142,2.111,2.079,2.046,2.011,1.975,1.938,1.9,1.86,1.82,1.778,1.735,1.69,1.645,1.598,1.55,1.501,1.451,1.399,1.347,1.294,1.239,1.184,1.129,1.072,1.015,0.958,0.901,0.844,0.786,0.729,0.673,0.618,0.563,0.51,0.459,0.409,0.361,0.315,0.272,0.231,0.193,0.158,0.127,0.098,0.074,0.052,0.034,0.02,0.01,0.003,0,0,0.004,0.012,0.023,0.037,0.054,0.074,0.097,0.123,0.152,0.183,0.216,0.251,0.289,0.328,0.369,0.412,0.456,0.501,0.548,0.596,0.645,0.694,0.745,0.797,0.849,0.901,0.954,1.008,1.062,1.116,1.17,1.225,1.279,1.334,1.388,1.442,1.496,1.55,1.603,1.655,1.707,1.758,1.808,1.858,1.906,1.953,1.999,2.043,2.086,2.127,2.167,2.205,2.241,2.275,2.307,2.337,2.365,2.391,2.414,2.435,2.454,2.471,2.485,2.497,2.507,2.515,2.52,2.524,2.525,2.524,2.522
];

L = [L L(1)];
U = [U U(1)];

E = 2.1e11;                 % 钢铁的弹性模量
rho_steel = 8456;           % 钢铁的密度
rho_oil = 864;              % 油的密度
rho_water = 1000;           % 水的密度
a = sqrt(E / rho_steel);    % 钢铁中声速
c = pi * a * gamma / (2.0 * sum(x));       	% 阻尼系数
grav = 9.81;                % 重力加速度

%% 计算
rho_fluid = rho_water*p_water + rho_oil*(1-p_water);    % 油井中流体(油水混合物)的密度
A = pi .* (R.^2);       % 各节截面积
T = 60 / chongci;       % 周期, 一个冲程的时间
w = 2*pi*chongci / 60;  % 圆频率

Wr = sum(A.*x) * (rho_steel - rho_fluid) * grav;        % 抽油杆各节总重

m = length(x);          % 抽油杆的节数
K = length(L);          % 采样点数量
nF = 10;                % 算几阶Fourier级数
t = linspace(0, T, K);  % 时间

%% Step 1
D = L - Wr;             % 动态载荷函数
sigma0(1)   = trapz(t, D.*cos(w*t*0))*w/pi;
tau0(1)     = trapz(t, D.*sin(w*t*0))*w/pi;
nu0(1)      = trapz(t, U.*cos(w*t*0))*w/pi;
delta0(1)   = trapz(t, U.*sin(w*t*0))*w/pi;
for n = 1:nF
    sigma(1, n) = trapz(t, D.*cos(w*t*n))*w/pi;
    tau(1, n)   = trapz(t, D.*sin(w*t*n))*w/pi;
    nu(1, n)    = trapz(t, U.*cos(w*t*n))*w/pi;
    delta(1, n) = trapz(t, U.*sin(w*t*n))*w/pi;
end

%% Step 2
for n = 1:nF
    alpha(n)    = n*w/(a*sqrt(2)) * sqrt( 1+sqrt(1+(c/(n*w)).^2) );
    beta(n)     = n*w/(a*sqrt(2)) * sqrt( -1+sqrt(1+(c/(n*w)).^2) );
    kappa(1, n) = (sigma(1,n)*alpha(n) + tau(1,n)*beta(n)) /    ...
                    ( E*A(1)*((alpha(n))^2+(beta(n))^2) );
    mu(1,n)       = (sigma(1,n)*beta(n) - tau(1,n)*alpha(n)) /    ...
                    ( E*A(1)*((alpha(n))^2+(beta(n))^2) );
                
    O(1, n) = (kappa(1,n)*cosh(beta(n)*x(1))+delta(1,n)*sinh(beta(n)*x(1)))*sin(alpha(n)*x(1)) ...
             +(mu(1,n)*sinh(beta(n)*x(1))+nu(1,n)*cosh(beta(n)*x(1)))*cos(alpha(n)*x(1));
    P(1, n) = (kappa(1,n)*sinh(beta(n)*x(1))+delta(1,n)*cosh(beta(n)*x(1)))*cos(alpha(n)*x(1)) ...
             -(mu(1,n)*cosh(beta(n)*x(1))+nu(1,n)*sinh(beta(n)*x(1)))*sin(alpha(n)*x(1));
    DO(1,n) = ( tau(1,n)/(E*A(1))*sinh(beta(n)*x(1)) + (delta(1,n)*beta(n)-nu(1,n)*alpha(n))*cosh(beta(n)*x(1)) )*sin(alpha(n)*x(1)) ...
             +( sigma(1,n)/(E*A(1))*cosh(beta(n)*x(1)) + (nu(1,n)*beta(n)+delta(1,n)*alpha(n))*sinh(beta(n)*x(1)) )*cos(alpha(n)*x(1));
    DP(1,n) = ( tau(1,n)/(E*A(1))*cosh(beta(n)*x(1)) + (delta(1,n)*beta(n)-nu(1,n)*alpha(n))*sinh(beta(n)*x(1)) )*cos(alpha(n)*x(1)) ...
             -( sigma(1,n)/(E*A(1))*sinh(beta(n)*x(1)) + (nu(1,n)*beta(n)+delta(1,n)*alpha(n))*cosh(beta(n)*x(1)) )*sin(alpha(n)*x(1));
end

%% Step 3
for j = 1:(m-1)
    nu0(j+1)    = sigma0(j)*x(j) / (E*A(j)) + nu0(j);
    sigma0(j+1) = sigma0(j);
    for n = 1:nF
        nu(j+1, n)      = O(j, n);
        delta(j+1, n)   = P(j, n);
        sigma(j+1, n)   = E * A(j) * DO(j, n);
        tau(j+1, n)     = E * A(j) * DP(j, n);
        
        kappa(j+1, n)   = ( sigma(j+1,n)*alpha(n) + tau(j+1,n)*beta(n) )    ...
                         /( E * A(j+1) * ((alpha(n))^2+(beta(n))^2) );
        mu(j+1, n)      = ( sigma(j+1,n)*beta(n) - tau(j+1,n)*alpha(n) )    ...
                         /( E * A(j+1) * ((alpha(n))^2+(beta(n))^2) );
                     
        O(j+1, n) = (kappa(j+1,n)*cosh(beta(n)*x(j+1))+delta(j+1,n)*sinh(beta(n)*x(j+1)))*sin(alpha(n)*x(j+1)) ...
                   +(mu(j+1,n)*sinh(beta(n)*x(j+1))+nu(j+1,n)*cosh(beta(n)*x(j+1)))*cos(alpha(n)*x(j+1));
        P(j+1, n) = (kappa(j+1,n)*sinh(beta(n)*x(j+1))+delta(j+1,n)*cosh(beta(n)*x(j+1)))*cos(alpha(n)*x(j+1)) ...
                   -(mu(j+1,n)*cosh(beta(n)*x(j+1))+nu(j+1,n)*sinh(beta(n)*x(j+1)))*sin(alpha(n)*x(j+1));
        DO(j+1,n) = ( tau(j+1,n)/(E*A(j+1))*sinh(beta(n)*x(j+1)) + (delta(j+1,n)*beta(n)-nu(j+1,n)*alpha(n))*cosh(beta(n)*x(j+1)) )*sin(alpha(n)*x(j+1)) ...
                   +( sigma(j+1,n)/(E*A(j+1))*cosh(beta(n)*x(j+1)) + (nu(j+1,n)*beta(n)+delta(j+1,n)*alpha(n))*sinh(beta(n)*x(j+1)) )*cos(alpha(n)*x(j+1));
        DP(j+1,n) = ( tau(j+1,n)/(E*A(j+1))*cosh(beta(n)*x(j+1)) + (delta(j+1,n)*beta(n)-nu(j+1,n)*alpha(n))*sinh(beta(n)*x(j+1)) )*cos(alpha(n)*x(j+1)) ...
                   -( sigma(j+1,n)/(E*A(j+1))*sinh(beta(n)*x(j+1)) + (nu(j+1,n)*beta(n)+delta(j+1,n)*alpha(n))*cosh(beta(n)*x(j+1)) )*sin(alpha(n)*x(j+1));
    end
end

for i = 1:length(t)
    cos_vec = cos([1:nF]*w*t(i));
    sin_vec = sin([1:nF]*w*t(i));
    
    u_xm(i) = sigma0(m) / (2*E*A(m)) * x(m) + nu0(m)/2 + sum(O(m, :).*cos_vec) + sum(P(m,:).*sin_vec);
    F_xm(i) = E*A(m)*( sigma0(m)/(2*E*A(m)) +  sum(DO(m, :).*cos_vec) + sum(DP(m,:).*sin_vec) );
end

figure; hold on;
plot(U, L, 'r-' ,'linewidth',2);
plot(u_xm, F_xm, '--' ,'linewidth',2);
legend('悬点示功图', '泵示功图','Location','Best');
xlabel('位移 (m)'); ylabel('载荷 (N)');
title('一级杆油井示功图_分三段_不独立');
grid on;