数学建模写作培训

赛题 要（19，并31Gibbs46针对问题四第一问中首先分析了Gibbs波动方程建立的过程认为Gibbs模型忽略了重力的影响，在Gibbs模型的基础上加入了重力因素加以改进，得到72得到了阻尼系数的求解，并给出了迭代求解阻尼系数的算法和迭代流程（图9：Gibbss模型，阻尼系数，傅里叶系数，有效问题目前，开采原油广泛使用的是有杆抽油系统（1。电机旋转运为部分圆弧、电机匀速运动，悬点E下只挂光杆（光杆下不接其它杆，不抽油，通常用来调试设备41出悬点E1Gibbs1井参数和悬点示功数据，附件2是有三级杆的另一油井参数和悬点示功数据，利问题三：泵功图的应用（22效冲程并利用附件1和附件2的数据分别估算两口油井一天（24小时）的产（C12问题四：深入研究的问题（2请对Gibbs模型进行原理分析，发现它的不足。在合理的假设下，重新附件1附件2的数据重新进行计算对计算结果与问题二的计算结果进行比较，Gibbs

a

c

其中数学方法估计参数c，贡献是很大的。对此，请你进行研究，详细给出计算c的理论推导过程并尽可能求出c。如果需要题目之外的数据，请用字母表示之并给出计算c模型通用符号说10为游梁初始摆2ww3HH为O4EE56II为O7SE8JJ为J对时间的导9为为游梁摆动的角速为任意时刻游梁摆动的角加速问题分

问题一杆悬点运动题目要求根据附录4给出的四连杆各段尺寸，利用附录1的参数，求出悬点1的悬点位移数据进行比较。E（位移函1有杆抽油系统四连杆运动函数计游梁摆动方程的建1，ABE对应坐标原点（称为E的下死点E0；D（时，E（称为E的上死点D（2）E0（即一个冲程）[1,2]。1BDLODJOOKODrH为OI为O到坐标在OBD中，根据余弦定理可得L2b2J22bJcos(得J2J2cos

Irsin2HrcosJ

tg1Irsin 可知式中wt

Hrcoscos1

J2L2 Irsin2HrcosJ

tg1Irsin Hrcos悬点位移函

为式中0为游梁初始摆角，即t0时刻游梁与纵坐标方向的夹角，故这里EA SE 1

J2L2 SESAa2cos

IIrsin2Hrcos tg1Irsin 悬点速度函

Hrcos将式（2）两边同时对时间t求导数，可得游梁摆动的角速度为JcosJbJsinJ和为J和对时间的导数。由式（3）和式（4）

Jwrsin(wrcos(J将式（8）和式（9）带入式（7）Jcoswrsin()wrcos(

bJsin

wrJsinbsinbJsin

JsinKsin

其中 K

tg1I2I2H将式(12)带入式（11）中，得游梁摆动的角速度wrKsinbsinbJsin游梁摆动时，因前、后臂的角速度相等，故任意时刻悬点的速度VE为V

=aawrKsinbsin 悬点加速度函

bJsinw2rKcosbcos2bwrcos2cos

bJsin(+

游梁摆动时，因前、后臂的角速度相等，故任意时刻悬点的速度WEWAaw2rKcosbcos2bwrcos2cosWE bJsin(+对于一定型号的抽油机其尺寸ab,L,H,I,r是已知的给定曲柄旋转的速度w或转速n，应用以上便可计算各运动参数[3]有杆抽油系统四连杆运动模很大误差，为此，本题中r/L14，这里我们把B点绕游梁支点的弧线运动近20时，游梁与连杆的连接点B在B处，为距曲柄最远的位置，相应于悬点E的下死点；时，B在B处，为距曲柄最近的位置，相应于悬E的上死点。在这种简化情况下我们对应的有[34]HOD+BDLIOB

1b2J2L2 SEa2cos

V

awrKsinbsin bJsin w2rKcosbcos2bwrcos2cos bJsin(+bbrsin2LrrcosJ br

Lrrcos b2b2Lrtg1

Lw a4.315

b2.495

r0.95

L3.675

n1YOB，COD，OOM，MODOAa，OBbBDLODJOOKODr；H为O到坐标横轴的距离I为O到坐标纵轴的距离SE为悬点位移函数；VE为悬点速度函数；WE为数据比3在一个周期内，利用（19）求得悬点处的理论位移（具体数据见附1111234问题分

问题二：泵功图计GibbsGibbs间变化的位移和荷载函数进而通过转换得到泵处随时间变化的位移和荷载Gibbs（位移函数Gibbs模型——一级杆情况下泵功图的计

c

其中：u(x,t)为抽油杆离地面悬点x深处在t时刻的位移；a为应力波在抽杆中的速度；c为阻尼系数波动方程边界条

u(x,t)|x0au(L,t)u(L,t)

初始条u(x,t)|x0F(x,t)|x0L(t)Wr

其中u(t)为实测光杆悬点位移函数mL(t)为光杆悬点荷载函数knWrn连续性条(Fi,j)1(Fi,j(ui,j)1(ui,j悬点位移和荷载方程的

D(t及光杆位移函数u(t展开成傅氏级数：D(t)2

NN(ncosntnsin

U(t)

(ncosntnsin22根据附件1和2中给出的悬点处的离散的位移和荷载的数据用下面的方法计算4个傅氏系数n,n,n,n，然后将傅氏系数带入到上面的悬点位移以n为例来说明求解过程：实际中的数据是离散的，需要将D(t)和u(t)离散化，设在悬点处进行时采样，采样点的个数为k144。然后令t，2，dtd

,N,N考虑到2k

,,

D(2p)

，U(2p)

Dcos(2n0)Dcos(2n

Dcos(2n1

cos(2n(k1))Dcos(2nk) kkkk2kkkkk2k

2D0cos0Dcos(2n1

cos(2n(k1))

cos(2n

k

k

,ND0Dkcos0cos2n,N2nnk

k

,,N2nnk

sin(k

nUp2nUpk

k

nUp2nUpk

sin(k

,N,N其中：为曲柄转角；为曲柄角速度；TN为傅氏级数所取项数，本文取10；p为从0到k各点序号k为2周期内等分的份数，即离散k144,N,NMATLAB画出悬点位移和荷载的曲线。粘滞阻尼系数求 C2

2(B

Aln 1)B1(L/a)sin(L/a)cos(L/a)m

ddqq d q q m2B12lnmBm4

(m2 d

lnA4E/E/ii其中：dy为油管内径，单位：m，d为抽油杆的平均直径，单位：mAi为m2；lmA为为抽油杆的平均截面积，单位：m2；为液体粘度，在本题中即为地面原油粘度，单位：Pa.sLmE为是抽油杆材料的弹性模量，单位：Nm3kgm3特殊函数和系数求x截面的位移随时间的变化。求解步骤如下[5,7,8]：首先令U(xt)X(x)T(tx和tT

CT(t)Xa2T

a2T Xn上式每一侧仅含有一个独立变量，因而它等于一个常数，令其为2，带n（39）nT(t)CT(t)2a2T(t)nn X(x)2X(x)n方程组（40）中（1）式的周期解为T(t)eint，带入方程组第一个等式可方程（40）的特解n11(c式子中n11(c

nn

nnn11n11(c

当n000，方程组（40）

T(t)CT(t)

X(x) T(t)X(x)方程组（40）综合以上各式，可求出方程（40）

X(x)nsinnxncosnu(x,t)(x)

sinxcos

u(x,t)

xx2

xF(x,t)

U(x,F(x,t)EA

O'(x)cosn

Pn(x)sinnt

u(x,t)

xx2

F(x,t)EA2EAOn(x)cosntPn(x)sinnt

On(x)coshnxnsinhnx)sinnx(nsinhnxncoshnx)cosnP(x)(ksinhxcoshx)cosx(coshxsinhx)sin On(x)EAsinhnx(nnnn)coshnx)sinn ncoshx()sinhx)cos

n n

Pn(x)EAcoshnx(nnnn)sinhnx)cosn

nsinhx()coshx)sin

n n k nnnk nn 上面各式中:E为钢杆弹性模量；Arx截面处的截面积；为多级杆情况下泵的位移和荷载计傅里叶系数的求 ,Am；各杆长度分别 ,ln；各杆的直径分别为d1,d2 ,dn；傅里叶系数为jn,jn,jn,jn，其j表示与系数有关的抽油杆段。类似的还有O(x),P(x),Ox),Pxj j j j j（j级抽油杆xj表jxjj

j

j0

j10xj1EAj

j

jnjn

j1On(xj1j0

j O' j j1j1 j P' 多级杆情况下泵的位移和荷载函

j j1j1 j5.1情况，在多级杆情况下，第j级抽油杆下端的位移和荷载函数为[7,8,9]：u(x,t)jj j

xj

2

NNjOn(xj)cosntjPn(xj)sin

j

N

F(xj,t)EAj2EAjOn(xj)cosntjPn(xj)sinnt

jOn(xj)(jkncoshnxjjnsinhnxj)sinnx (sinhxcoshx)cos j n j n nP(x)(ksinhxcoshx)cosj j n j n n (jncoshnxjjnsinhnxj)sinnxO(x)jnsinhx()cosh

x)sinj

n jn jn nj n jncoshx()sinhx)cos

n jn jn n nP(x)

j

cosh

(

)sinh

x)cosj

n jn jn nj n jnsinhx()coshx)sink

jnnjn

n jn jn n nj

j

jn jn 上面各式中:EArx计算流本文在所给悬点位移和荷载数据的基础上，通过利用Gibbs模型，推导出求泵的位移和荷载的方法及详细运算过程。下面给出计算流程步骤及计算流程图[7]步骤1k份，即离散化的采样点数，本文取k144（34（374Dt)和Ut)的函数，即可获得悬点处位移和荷载的曲线。步骤2：计算粘滞阻尼系数C。利用（38）计算阻尼系数步骤3:计算特殊系数和函数系数包括n，n，jkn，n函数包括Qn(x) 计算knn，Qn(xQ(xP(xPn(x 56。jknjnjQn(xjjPn(xjjQn(xjjPn(xj4。F(xt)和位移函数U(xt)，进而可得出井下示功图。是否为最后否是否为最后否是44计算结2v---------------------37位移荷载注：具体数据见所提 .xls文5741位移荷载注：具体数据见所提 .xls文问题分

67问题三：泵功图的应用——泵内充气判本题是对泵功图的应用，本文选择对第2（泵内气体判断）进行研究。气体影响情况分这里我们换法再次确定泵功图的塞柱有效冲程我们知道有杆泵主要由泵筒,柱塞，游动阀和固定阀组成（在这里，我们用TV来表示游动阀，用SV来表示固定阀如下图7所示）。所以我们根据泵功图来确定阀开闭点的位置。同PS[10,11]7若有气体进入泵腔，柱塞下行时气体将受到进一步压缩，使TV延迟开启。设柱塞位于上死点，泵腔内气柱的高度为Hg，压力为Pi。假定柱塞下行时，气HHgHPP HgHgH

PHH1PP0，使TV开启，即P

H，并带入式(61)Hg

H1

P

PiH

11PH 0 1HH1时，TVTVPP0保持不变。当柱塞从下死角开始上行时，泵腔内的压力P减少，泵腔内的气体将PPiSV开启。降低）P的变化趋势是一致的。即泵腔内气体的存在将使柱塞在上、下冲程阶段缓慢地加载或卸载；SV、TV都将延迟开启。但是，只有下行程TV的延迟开启才影响泵效。TV开启后，柱塞下有效冲程的计SVTVP随S的增大而减小，FP

K0SV关闭后，P将随S的增大而增大，FP则随之减小。泵示功图这一曲K0。下行程TVSVP随STV开启后PP0为一常数FP表示不变。泵示功图上这段曲线为一水平线，K0。在TV关闭后，PS的减小而减小，FPK0。K发生突变。因此，通过第二题，我们可以得到泵示功图的一组离散数据点PiXi0,Yi0（i n）iiXXi0XLii

Yi0

XRX YmXRXL表示泵示功图左右两端点的横坐标；Ym,Yn分别表示泵的最大Pi点前后各取k（k2k1Pi点前的斜Ki

iX2Xi

Xi

j

j2k1,Yj

jjikj

2kii

jik

X

j

22kijjKi后，结合泵功图的几何形状和几何意义，可按下述依ijjSV开启点：在上行程（Xi1XiPi点前的k个点中至少有k1点的斜率KPi点后的k个点中至少有k1个点的斜率K;则认为Pi点SVPiP1。SVPikk1KPi点后的k个点中，至少有k1KPiTV开启点：在下行程（Xi1XiPi点前的k个点中至少有k1个点的斜率KPi点后的k个点中至少有k1个点的斜率K;则认为Pi点就是TVPiP3'。TVPi前的k个点中至少有k1KPi点后的k个点中，至少有k1KPi点就是TV关闭PiP4。SV和TVsmins1,s2

s1|P2'P1'||X20'X10'

s2|P4'P3'||X40'X30'泵内充气判步骤1：首先通过求解吉布斯方程，可以求解得到一组离散的数据点PiXi0,Yi0),(i n步骤2：在Pi前后各取k2k1并求出Pi点处的斜率Ki

iX2Xi

Xi

j

2k1,Yi

jik

ijjijj

jik

X

j

22kijj步骤3：求得各点处的斜率Ki4ijjSVP1，并记录此时的0.3P1前的k有k1KP1后的k个点至少有k1KSVP2，此时的0.3P2前的k个点至少有k1个KP2后k个点并没有k1K。找到TVP3P3前的k个点并没有k1KP3后的k个点至少有k1K找到TV的关闭点记为P4P4kk1个点的斜率KP4点后的k个点中，至少有k1K计算结2）的有效冲程、泵内是否有气的判定结果为：5有效冲程（米有问题四：深入研究的问问题分题目要求对Gibbs模型进行原理分析，发现它的不足。在合理的假设下，重的详细计算过程。第二问要求计算的理论推导过程并尽可能求出首先分析Gibbs波动方程建立的过程得到Gibbs模型忽略了重力加速度的影响，因此我们在原有的基础上对Gibbs模型加以改进，增加抽油杆的重量对模Gibbs首先分析Gibbs波动方程建立的过得到了Gibbs模型忽略了重力加速度的影响，因此我们在原有的基础上对Gibbs模型加以改进，增加抽油杆的重量对Gibbs模型改Gibbs模型分设符号x是悬点移动前油杆的位置，取正下方向下，取u(x,t)作为位移

fxfxxfamfd作为单元体单位长度上粘性阻力尼，其方向与速度方向相反fw取为单元体重力，它们的单位都是牛顿，除此之Er代表抽油杆的材料弹性模量，单位是kg3；vem粘滞阻力系数，单位是kgm.s根据受力分析可知，作用在单元体上的要有一下几个fEA(u rrx EA(u r

farAr(2tfv(u rfA(u rr根据单元体轴向力平衡条件Fx0 ErAr[((x)xx)(x)x]rAr(t2)xvetxrAr即

2urrx2

rAr

vuAe r

2 t2

x2c

am改进Gibbs模

a c为阻尼系数s1c

rGibbsGibbs设浸入在油管中第J2uJ(x,

2uJ

uJ(x,t)

0；a为应力波在抽油杆中 速度；g为重力加速度0EAJJ段油杆的截XJJ段油杆的长度。方程（1）uJ(x,t)

QJ(x)cosntPJ(x)sinntgfx2BJx

QJ(x)(kJcoshxJsinhx)sinx(JsinhxyJcoshx)cos PJ(x)(kJsinhxyJcoshx)cosx(Jcosh n

11(c/ a

1(c/ akJ

nnn EA(22 nn

EA

22 BJ

AJJJJJ

BJ1FJJJAJJJ

aJ1GJbJ1H bJ1GJaJ GJcos(x)sinh(x)sin(x)cosh( HJcos(x)sinh(x)sin(x)cosh( a

AJ1eJ

BJ1

J

aJ1

JbJ1h J

J

b

BJ1eJ

AJ1

J

bj1

JaJ1h J

J

sin(x)cosh(x)cos(x)sinh(eJ f

nsin(nx)cosh(nx)ncos(nx)sinh(n J J

)2

g

x)cosh( hJsin(x)sinh( 根据式（71）J

uJ(x,t)计算结

F(x,t)

27（28Gibbs7阻尼系数求抽油杆中长度为dx的单元体的粘滞阻尼力为FCudxc

u F

r其中c为粘滞阻尼系数为抽油杆的密度；Ar为第r级抽油杆的横截面理接下来的几个中的产液量参量也当作已知量。 S2S

T1000

利用式上式可推导阻尼系数的计算如下225.08(PP)Tmc h m

S

AsgtN

AsgtNmN1sP1.135107QL代表水的密度，Kgm3；L代表有效提升高度，m pn Q代表实际产液量，m3d；L

1000PtPc；P为套压，Pa；Pct ct

L 1.67105 n； 为泵功图的面积 Nm；

1440ASPN；A为泵（油管）m2

=S- S mSm0

n pLmA为泵（油管）m2Kgm3 bbF缩b(1F)b wrhP代表光杆功率，W；代表抽油杆密度，Kgm3；T代表抽油机的一个冲程时sg代表重力加速度，9.8m/s2P代表水利功率，W。rh阻尼系数迭代求解算法步骤为步骤2计算阻尼系数c步骤3：Gibbs波动方程，求解泵示功图；步骤4：根据泵示功图计算有效冲程SPe6步骤6：根据泵示功图求泵功率Ppump7：P

Ppump，是则该c即为所要求的数值；否令 转

Ph6。其中1和2力系数c。算法流程图如下所示： NYNY9参考文王鸿勋,张琪《采油工业原理[M].石油行业,1985年4月,第2章万仁溥.《采油工程手册》[M].石油行业,2000年8月,第5章,第节:365-王常斌,陈涛平,郑俊德.游梁式抽油机运动参数的精确解[J].向刚.稠油井杆中管掺热流体闭式循环举升工艺设计[D].大,2011:38-Gibbs.S.G,MethodofDeterminingSuckerRodPumpPerformance.UnitedStatesPatentOffice,Sep.1967. 张琪, .抽油井计算机诊断技术及其应用[J].华东石油学院学报2:144-周继德.用示功图拉计算抽油机井产量和泵筒漏失量[J].油田地面工程,1986,(2):60-62.朱君.有杆抽油系统井下工况诊断方法研究[D].大庆石油学院博士 puterDiagnosisofDownholeConditioninSuckerRodPum 附录2u1Columns1through - - 11through21through31through41through 51throughColumns61through Columns71through Column-fColumns1through Columns11through Columns21through Columns31through Columns41through Columns51through 61through Columns71through Column3：ylabel('m');holdonholdoffforforforforafunctionu=beng(x,t)kesi=[89.3227-7.4966 -0.2615- -0.00382.66E-04- -0.0299-0.0248 -0.00122.64E-05-7.93E- 2.42E-051.12E-%v=[71.5291 -0.3797 -0.7052 -0.5084 - -%deta=[-23.8556 -2.3054 -0.5375 -0.2668- -0.9275 -0.2159 -0.1154 -0.2669 deta=[- -3.0422-0.2337- forn=1:1:10