2012年全国数学建模竞赛优秀选

（由填写题 要180360°曲柄转问题二：结合参考文献建立了泵功图的Gibbs模型，详细推导和讨论了悬点示件的定义和物理意义,并用编程语言对模型进行了实现,绘制了完整的程c的计算进行了详细的理论推导，获得了最终的计算表达式。此外,利一 问题重 二 问题假 三 符号说 四 问题一：光杆悬点运动规 问题分 模型建立及推 悬点运动参数方 悬点运动速 悬点运动加速 模型初始条 阈值范围及初始条 模型验证及误差分 结 五 问题二：泵功图计 Gibbs模型的基本思 Gibbs波动方程的推 Gibbs波动方程的求 一级杆Gibbs波动方程的求 多级杆Gibbs波动方程的求 模型详细实现过 原始数据的处 边界条 初始条 计算步 程序设计流程 模型的验 一级杆油井的模型验 三级杆油井的模型验 小 六 问题三：泵功图的应 背景分析与假 问题一有效冲程法模 模型的建立与推 模型的验 产量计 问题一有效功率求解法模 模型的建立与推 模型的验 产量计 问题二模 模型的推 模型的建 6.4.5模型使 七 问题四：Gibbs模型深入研 问题 基于摩擦功的阻尼系数模型建 模型求 八 模型评 九 参考文 一、题Gibbs模型，给出由悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；根据附件1、3油井参数和2个不同的由泵功图估计油井产量的模型，至少一个要利用“有效对Gibbs模型进行原理分析，发现它的不足。在合理的假设下，重新建立抽2u2 Gibbst2ax2cta已知常数，c称为阻尼系数，鉴于大多数的阻尼系数是作了诸多假设后推出数c，详细给出计算c的理论推导过程并尽可能求出c。二、题假三、号AO长 曲柄半支点与曲柄滑块D 支点O与曲柄圆心O’间距离 OO’水平 任意时刻J和-Y方向夹角 四、问题一：光杆悬点运动规律问题分析t模型建立及推导4.1运动，在ΔOBD中，根据余弦定理L2b2J22bJcos[(其 Jtg1 Ir HGr

(4-2-(4-2-(4-2-式中，t。由式（4-2-1） 1b2J2L2cos (4-2-4-2-4J和都是只为的函数，因此也只为悬点运动参数方程位移对于任意时刻，后臂的角位移 (4-2-式中，0为后臂的初始摆角，即t=0时刻，后臂和Y方向夹角，可令t=0（4-2-求得。后臂绕指点OSA[0 (4-2-一常数，因此SA只是的函数，给定值，便可以计算出悬点位移SA。悬点运动速度将式（4-2-4）两边同时对时间t求导数，得到后臂摆动的角速度'Jcos()J'bJsin(

(4-2-（'rcos(J

(4-2-(4-2-(4-2- 4-2-9)代入到（4-2-7）中，整bJsin(

Jsin()Ksin(0 K

0tg H（4-2-11）（4-2-10）'r[Ksin(0)bsin(bJsin(

,故任意时刻悬点的速度VAVa'ar[Ksin(0)bsin( bJsin(

悬点运动加速度4-2-12式两边同时对时间t2r[Kcos()bcos()]2br'cos()'2cos('' bJsin(2r[Kcos()bcos()]2br'cos()'2cos(Wa'' bJsin(式（4-2-5）（4-2-12）(4-2-14)反映了后臂摆角的变化规律,（4-2-6）,(4-2-13)和式(4-2-15)是悬点运动的位移、速度和加速度计算。对于一定型号的抽油机,其几何尺寸（a,b,L,H,G,r）已知的，给定曲柄旋转的角速度或转速,应用以上便可计算各运动参数。模型初始条件1.始时刻0，代入计算得到初始角度01190。因此为(4-2-6)可以表S[1200 阈值范围及初始条件O'点与O点相对位根据题设条件建立的抽油杆三个时刻运动模型（OB水平线，曲柄滑块在顶点D，后臂位于BD’曲柄滑块处于最底端，此刻后壁与连杆间张角最大（<180），O点到滑竿尾端的距离等于OD’’D’’D点。为保证连杆的正常运动，D’’bLOD''bLOD为保证抽BM2MD2(BMOB)2(MD2r)2(bBM=3671.7mm,ML=155.4611mm。所以为保证连杆系统的正常工I2614.4611,H4621.7。4.2O'B部分资料中提到O'BOO'2OM2MD2=此时OD=OD''2OM2MD''26107.83mm2此时ODOD''2OB2BD''22OBBD''cosOBD''OB2BD''261702所以，得到BM1664.54715mm，即I （b）O'位于B模型验证及误差分析用实现推导的算法，建立有杆抽油系统模型，图4.4表示悬点运动示功图4.4（a）是对抽油杆选点出荷载变化与位移变化规律及相互关系的直观描述。用计算得到模型的悬点示功图。A是悬点运动的下点，C是悬点运动的上死点，从ABC点的过程，是悬点的上行；CDA'AB，油管卸载而缩短，抽油杆的静态形变完毕，柱塞开始相对于泵向上移动，位移增大。由C到D过程中，抽油杆卸载而缩短，油管加载拉长，悬点位移变化较小；由D到A'过程中，油管与抽油杆的静态变形完毕，柱塞开始相对于泵筒向下移动。1803604.4(c)，表示有杆抽油机的一个工作冲程中，下冲 图 (a)悬点位移曲线模拟结果和附件1实际测试曲线拟结180360°曲柄转角，通过计算五、问题二：泵功图计算Gibbs模型的基本Gibbs给出了悬点示功图转化为示功图的模型。如果有杆抽油系统中的抽油杆看作一段传输线，并把抽油杆上端点看作，把抽油杆的下端杆以一定波速(波速的大小由抽油杆的密度和弹性模量决定)传输到抽油杆上端据处理，即获得位移时间图和荷载时间图，为了实现这个过程，Gibbs60年Gibbs波动方程的推由牛顿第二定律可以知道，F=ma，对于抽油杆而言，它受到的要有抽律可以得到张力计算：𝐹1(x,t)=𝐸𝐴∙ ,𝐹2(x,t)=𝐸𝐴∙ 其中：𝐸6可知𝐸2.11011(pa);A为抽油杆的截面积;𝑢(𝑥𝑡)𝑥1为抽油杆上端点坐标，𝑥2为下端点坐标，且有𝑥2𝑥1L，L为抽油杆长 可以知道 𝜕 + 𝜕 +L 现取前两项作为数值近似计算，可以得𝐹(x,t)=𝐸𝐴∙ ≈𝐸𝐴∙ +L∙ 3𝐹(x,t)=𝑐∙𝑚𝜕𝑢(𝑥,𝑡)=𝑐𝜌𝐴𝐿 3

𝐹1(x,t)−𝐹2(x,t)−𝐹3(x,t)=𝐸𝐴∙L

𝜕2𝑢(𝑥,|

𝑥

𝜕𝑢(𝑥,=𝑚𝑎=𝜌𝐴𝐿2

（5-2-对上式进行运算简化之后即可以得到Gibbs𝜕2𝑢(𝑥,𝑡)=𝑎2𝜕2𝑢(𝑥,𝑡)−𝑐𝜕𝑢(𝑥,𝑡) 第一类边界条件:𝑢(x=0𝑡)=第二类边界条件

𝜕𝑢(x=0,𝑡)

D(𝑤t)

{初始条件𝑢(x0𝑡0U(0)𝑢(x0𝑡0其中：𝑎=√𝐸为应力波在抽油杆 的速度，单位为(m/s)，由附录6知抽油杆的密度为𝜌=8456(kg/m3)，所以𝑎=√2.1×10c为阻尼系数𝑢(𝑥𝑡)由(5)可以知道，Gibbs模型的本质是利用牛顿第二定律和传输线理论构Gibbs波动方程的一级杆Gibbs波动方程的求解为了进一步求解Gibbs波动方程，可以用复变量z(𝑥𝑡)代替方程（5-2-6）中的𝑢(𝑥𝑡)，则有：𝜕2z(𝑥,𝑡)=𝑎2𝜕2𝑧(𝑥,𝑡)−𝑐 𝑧(𝑥,𝑡)=𝑋(𝑥)𝑇(𝑡) 关于x，𝑇(𝑡)关于t𝑋′′(𝑥)=𝑇′′(𝑡)+ 𝑋′′(𝑥)=−𝜆2=𝑇′′(𝑡)+ 𝑋′′(𝑥)+𝜆2𝑋(𝑥)= 𝑇′′(𝑡)+𝑐𝑇′(𝑡)+𝜆2𝑎2𝑇(𝑡)=0 为了确定参数𝜆的值，可以对𝑇(𝑡)关于时间t进行周期傅里叶级数展开(𝑇𝑛(𝑡)=𝑒𝑖𝑛𝑤𝑡)，然后代入到（5-3-6）式中确定参数𝑒𝑖𝑛𝑤𝑡(𝜆2𝑎2−𝑛2𝑤2+𝑖𝑐𝑛𝑤)=由于𝑒𝑖𝑛𝑤𝑡≠0，而上式要成立，则必有𝜆2𝑎2−𝑛2𝑤2+𝑖𝑐𝑛𝑤=0定𝜆n的值𝜆=−α+

1 1+ −1 1+ 𝑛𝑤)+𝑖 个基本解是有正弦函数和余弦函数线性组合迭加形成的，结合上式的求得的cos(𝜆n)5式的常量解可以表示为T(𝑡):ξ,𝑋(𝑥):𝜂𝑥,𝜁，因而复位移函数可以表示为：∞𝑧(𝑥,𝑡)=𝑋(𝑥)𝑇(𝑡)=ξ𝜂𝑥+ξ𝜁+∑(Φ𝑛sin(𝜆nx)+Θ𝑛其中：Φ𝑛=−𝑘𝑛−𝑖𝜇𝑛和Θ𝑛=𝑣𝑛−𝑖𝛿𝑛都是复常量。首先，利用（5-2-6）式中的第二类边界条件可以有𝜕𝑢(x=0,𝑡)=D(𝑤t)= Re(𝜕𝑧(x=0,𝑡))，那么由z(x,t)

nncosn

nnsinD(t)EA[ (cosnxnnsinnx)eint nn

nEA(nnnn),nnEA(nnnn),n0

nn EA(22),n EA(2 nnnnEA(2 这时，为了满足第二类边界条件参数n和n需满足如下两个等式 D(t)cosntdt,n D(t)sinntdt,n U(t)cosntd,tn2n2 U(t)sinntd,tnn2其中，n和n

给出，且0

Re[z(0,t)]U(t)。

(ncosntnsinu(x,t)0x

n[On(x)cosntPn(x)sin

On(x)(ncoshnxnsinhnx)sinnx(nsinhnxncoshnx)cosnPn(x)(nsinhnxncoshnx)cosnx(ncoshnxnsinhnx)sinn) F(x,t) (On(x)cosntPn(x)sin

O'(x)[nsinhx()coshx]sinx n n [ncosh

x()sinhx]cos

n n P'(x)[ncoshx()sinhx]cosx n n [nsinh x n )

x]sinx这样，我们就完成了单级杆Gibbs方程的求解多级杆Gibbs波动方程的求L1、L2和fr1及fr2分别为第一级和第二级抽油杆柱的长度和截面积。用移函数u(L1,t)和荷载函数F(L1,t) u(L,t)1

L 21 21

[1On(L1)cosnt1Pn(x)sin

1 F(L1,t)Efr1[ (1On(L1)cosnt1Pn(L1)sin

10L O(L P(L22

11

2n1 2n1

1Pn20

2

2

2 u(x,t) x220

[O

)cosnt

2 2Efr2

22 F(x2,t)Efr2[ (2On(x2)cosnt2Pn(x2)sin

r 式中，2On(x2)及2Pn(x2)的计算与单级杆时的相同，如需要采用新的模型详细实现过程原始数据的处理【2GibbsGibbs【3Gibbs模型的初始条件，为了便于傅里叶系数的求边界条件Gibbs一维粘滞阻尼方程的推导可以知道，Gibbs模型主要包含两类边界第一类雷边界条件:𝑢(x=0,𝑡)=第二类纽曼边界条件

𝜕𝑢(x=0,

𝐸6可知𝐸2.11011(pa)；A为抽油杆的截面积。初始条件𝑢(x=0,𝑡=0)= 𝑢(x≠0,𝑡=0)=计算步骤本文对于Gibbs 2KnK

Dp

),n0,1,2,...,

2KnKUU

Dp

),n1,2,...,U2U

2n

K),n0,1,2,...,K

2

2n

),n1,2,...,KK KK1

1 a11ca11cnwKa1a11cnw【3nn

,n

,n1,2,...,nnnnnEA(2

nnnnn,n1,2,...,EA(2 【4】利用以上两步获得到的系数计算任意一点处位移函数和荷载函数的系数O(x)P(x)Ox和P'(x)； 【5】计算位移函数u(xtF(x,t) EAO'(x EAP'(x 0 j1 j j0xj O ), Px 0j j 程序设计流程图计算系数计算系数，，和系数𝜎𝑛𝜏𝑛𝑣𝑛,对傅里叶待定系数𝜎𝑛,𝜏𝑛,𝑣𝑛,𝛿𝑛模型的验证一级杆油井的模型验证4x悬点悬点泵功

一级654荷载荷载210---- 5.1泵功图是对抽油杆柱最下端荷载变化与位移变化规律及相互关系的直观描（51图有以下几方面的不同：F(x,t)EA[0(O'(x)cosntP'(x)sinn

，其中n=仿真求解过程中，n取值不可能取到，故用一个较大的数来代替，因三级杆油井的模型验证5.2三级杆悬点示功图和泵功图小通过上述分析，比较泵功图的拓扑相似度、几何形状，Gibbs模能比较准确六、问题三：泵功图的应用背景分析与假设第二问的模型一：有效冲程法模型模型的建立与推导n节，每节长分别为l1l2ln m1Al (Sm1l 0im,0m

pi m

Sl其中Api理论上应为对应各节的油管横截面面积V式中，V：实际产液量，单位为m3次Vth：理论产液量，单位为m3次sfs为柱塞冲程系数，是柱塞冲程Sp与悬点冲程Sr之比，即：sSpSr，f为泵筒充满系数，是柱塞有效冲程Spe与柱塞冲程SpfSpeSp lSpe计算得到的井下排液体积为V，漏失体积为V，则漏失系数V Vpe，其对应的地面体积为Ve，则体积系数vVeVpe。 V(ApiliApm(Sli)Vp) Vp：抽油泵一个冲次的漏失量，单位为m3 ：泵排出压力条件下油管内混合物的体积系数。计算柱塞有效冲程S计算抽油泵的漏失量Qp计算泵排出压力条件下油管内混合物的体积系数 柱塞有效冲程据做归一化处理，归一化后的荷载比记为Rf，位移比记为Rs。曲线转化为在区间[0,2]上的单值曲线Fs。泵功图曲率计算模型（五点差值泵功图曲线上任意一点PiKi根据与其相邻的五个数据点：Pi2(si2,fi2)Pi1(si1,fi1Pi(si,fiPi1(si1,fi1Pi2(si2,fi2之间的621Pi点的曲率为：KK i式中，i为直线Pi2Pi到直线PiPi2的有向旋转角度，计算为(fifi2)(sisi2)(fifi2)(sisi2)iarctan 弧长li的近似计算如下

1(sisi2

si2 liPi2Pi1Pi1PiPiPi1Pi处的曲率变化量KiPiPi1两点处的曲率变KiKi1KiKi2Ki1KiKi1Ki2/Sfoc和游动阀开闭点的位置差Smoc，则柱塞有效冲程Spe等于其中较小值归一化后的泵功图由区间S[0,1]上的闭合曲线转化成在区间S[0,2上的单值曲线F(S)按照等间隔S在区间S(0,1)，F[0.8,1.0]上依次找出曲率变化最大的两个点Eg1sg1,fg1Eg2sg2fg2，则可得比SRgsg2sg1S1,2)F[0,0.2]上依次找出曲率变化最大的两个点Ey1sy1,fy1Ey2sy2fy2，则可得比例下SRysy2sy1位移比Rs即可得到实际有效冲程模型的验证一级杆有效冲程计算步骤：移比记为Rs。原始泵功图的荷载归一1

原始泵功图的位移归一1 归一化荷载归一化荷载比归一化位移比 0 采样

0 采样6.2.2的单值曲线F(S如图6.2.3。归一归一化泵功图的区间归一归一化荷载比0

6.2.36.2.4均匀采样后的单值曲线单值单值曲线F(S)的均匀采归一归一化荷载比0

6.2.4一级杆均匀间隔采样F（s）6.2.4Eg1(sg1,fg1Eg2(sg2,fg2，求SRgsg2sg10.812 y 点Ey1(sy1,fy1)和Ey2(sy2,fy2)，则可得比例下位移 y

0.816SRg与SRy0.812Rs2.582.0951曲率曲率变化量0 平均曲率变化量平均曲率变化量Y:0.2Y:

Y:

Y:FF归一归一化荷载 6.2.5三级杆有效冲程计算步骤：1归一归一化荷载比

原始泵功图的荷

1归一归一化位移比

原始泵功图的位0

采样点6.2.6将归一化后的泵功图由区间S[0,1]上的闭合曲线转化成在区间S[0,2]1归一归一化荷载比0

6.2.7

均匀采样后的单值曲线1归一归一化荷载比0 扩展区

6.2.8三级杆均匀间隔采样F（s）Eg1(sg1,fg1Eg2(sg2,fg2，求得SRgsg2sg10.668。Ey1(sy1,fy1Ey2(sy2,fy2SRysy2sy10.606比较SRg与SRy0.606，对其乘以归一化位Rs5.1933即可得到实际有效冲程3.147米。压缩压缩比:曲率曲率变化量0 扩展区平均曲率变化平均曲率变化量 Y:

压缩比: Y:

Y:

X:1Y:

归一化荷载归一化荷载比0 扩展区

6.2.9产量计算利用模型求产油量QQNs第二问的模型二：有效功率求解法模型模型的建立与推导马力，用符号PRHP表示，且由下式计算：NsttPRHP

60

NsFS(t)：表示悬点位移值，单位为mNs：表示抽油机的冲次，单位为min1PRHP代替，这样光杆功率就可用泵功图所圈闭的面积与光杆运动周期之比。擦功，又称水马力，用符号HPh表示，由下式计算HP gQH Qm3He：表示有效举升高度，单位为m力所做的功，又称摩擦马力，用符号HPf来表示。抽油杆柱与油管的摩擦力：在直井中通常不超过抽油杆重量的1.5%（1（2）向下，故增加悬点荷载。下冲程中作用在悬点上的摩擦荷载是由（1（2（3）抽油杆柱与液柱之间的摩擦力大值可由下面的近似来确定： m2 f 9.81103(m21)Lm(m21) 式中，Ff1——抽油杆柱与液柱之间的摩擦力，单位为L——抽油杆柱长度，单位为——井内液mmdtdrdrVma——抽油杆柱最大下行速度Vmax可以按悬点最大运动速度来计算，计算时采用下面的近似（把悬点： S 式中——曲柄角2nS由上述可以看出，决定Ff1的主要原因是井内液体的粘度及抽油杆液柱与油管之间的摩擦力 Ff液体通过游动阀时产生的摩擦1V f2Vh f2

ff2fVfgpf——活塞截面积，单位为cm2pff——阀门孔截面积，单位为cm2 ——由实验确定的阀流量系数，对于常用 阀可由下图查出

df 式中，NRe——雷诺数dfVf——液体运动粘度；单位为cm2s

V(t)Ssin VS 将V代入式h

f2V

f f

f2

ph f2p7292f f3F f f 7290f

f F——液流通过游动阀所产生的活塞下行阻力，单位为1——液

(Ff1Fu1F关系式就可得到油井的日产液量Q模型的验证产量计算12问题二模型模型的推导的面积比例，以充气程度达到kj(0kj1)为例：kSA其中j则

D为理想泵功图面

A为充气泵功图面积SS AS'S'S化率之比S，由泵功图易知S1，将此值作为一个评价标准，若有当前泵模型的建

ttNsF(t)S(t)dt(F(T/2)F(0))(S(T/2) 60SA其中，为充气百分Kj接下来，我们考虑实际情况下，没有额外损耗的泵功图。这可由泵效 给是由于影响泵效的因素有很多，矿场实践中平均泵效往往低于0.7。因此为了DDDD

S' S

k

Ak或或S

S' k

Ak比S，结合前面的讨论，可以简单认为当S'1时，泵内有气体的结论成立。6.4.5模型使1归一化荷载比归一化荷载比0 6.4.1均匀采样后的单值曲线单值单值曲线F(S)的均匀采归一归一化荷载比0

6.4.2s区间[0,1]s区间[1,2]曲线包含的面积即为归一阀值0.7阀值0.6阀值阀值归一化荷载比1归一化荷载比10016.4.31归一归一化荷载比0

6.4.4题四：Gibbs模型深入研问题基于摩擦功的阻尼系数模型建立首先，长度为dxF

udxc

u F r

（7-1-则长度为 的抽油杆柱在一个循环中总的粘滞摩擦功可以表示为2L LTu 2Wo0Fdtdtdxcrfr00tdtdx（7-1- c rTu

（7-1-rfr00t LTu LT W2 m 1(B11) dtdx (B11) 00t 0 xL （7-1- 1v gp g （7-1-r rr

x 1v 1

p （7-1-r rrr g

xdv0Apg

rvA

x

（7-1-

vAr2Clnr 利用边界条件 u（抽油杆柱的运动速度）和 0可得C和C