机械系统动力学报告游梁式抽油系统曲柄运动规律的动态仿真

1、机械系统动力学项目指导书 课程三级项目报告游梁式抽油系统曲柄运动规律的动态仿真 课程名称：机械系统动力学 任课教师：XXX 班 级：XXX 团队成员：XXX 2016年6月44目 录一 项目背景0（1）游梁式抽油机简介0（2）我国抽油机现状1二 力学与数学模型1（1）抽油机运动规律建模2（2）曲柄运动规律建模7三 数学仿真方法12（1）数学工具简介12（2）龙格库塔法简介12（4）部分代码13四 系统动力特性仿真分析16五 收获与感想22六 参考文献22七 附录（程序代码）23一 项目背景（1）游梁式抽油机简介游梁式抽油机是油田目前主要使用的抽油机类型之一，主要由了驴头游梁连杆曲柄机构、减速箱

2、、动力设备和辅助装备等四大部分组成。工作时，电动机的传动经变速箱、曲柄连杆机构变成驴头的上下运动，驴头经光杆、抽油杆带动井下抽油泵的柱塞作上下运动，从而不断把井中的原油抽出井筒。主要特点：1、整机结构合理、工作平稳、噪音小、操作维护方便；2、游梁选用箱式或工字钢结构，强度高、刚性好、承载能力大；3、减速器采用人字型渐开线或双圆弧齿形齿轮，加工精度高、承载能力强，使用寿命长；4、驴头可采用上翻、上挂或侧转三种形式之一；5、刹车采用外抱式结构，配有保险装置，操作灵活、制动迅速、安全可靠；（2）我国抽油机现状目前，我国抽油机主要制造厂有十几家，产品主要以游梁式抽油机为主，约占抽油机总数的98%至99

3、%,有30多种规格，并已形成了系列，前置式抽油机，异相曲柄平衡抽油机，前置气平衡抽油机，配有CJT型节能拖动装置的常规型抽油机，和用窄V带传动的常规型抽油机等均与在各个油田推广应用，并取得了显著的经济效益。常冲程、低冲次的无梁式抽油机的研制也取得了也一些进展，如胜利油田设计并并与有关厂家协作生产的链条式常3冲程抽油机，已有近千台在各油田投入使用，在低冲程抽油机和抽稠油方面已初见成效。此外，桁架结构的滑轮组增拘束抽油机，滚筒式长冲程抽油机进入了实用阶段；次轮增距式长冲程抽油机的研制工作也取得了新的进展。重量轻、成本低、便于调速和调整冲程的液压抽油机，经过几年的研制和工业性试采油也积累了一定的经验

4、，其它形式的新颖抽油机加代传动游梁式抽油机，新型遥杆抽油机，大轮式游梁抽油机，留连干游梁式抽油机也正处于开发和研制过程中。然而，游梁式抽油机的缺点是不容易实观长冲程、低冲次的要求，因而不能满足稠油井，深油井和含气井采油作业的需求。同时长冲程、低冲次的游梁式抽油机尚待完善，（如油田正在使用的链条式抽油机还存在链条寿命短，换向冲击载荷大河钢丝绳易断，导轨刚度不足容易变形等问题）而且品格还很少，不能适用于当前石油工业的发展，液压抽油机至今还处于研制阶段。二 力学与数学模型为了对异相曲柄平衡抽油机的优化设计，以放映抽油机动力性能的最大扭矩因素为最小目标函数优化抽油机的机构尺寸，将整个抽油机装置简化为以

5、曲柄为等效构件的运动情况，我们做了以下假设条件：1、 电网供电电压与供电频率是常数；2、 不考虑电动机转子到减速箱曲柄轴各传动副的间隙与传动件的弹性变形；3、 在分析抽油机主体机构的运动规律时，不考虑各传动副的间隙与传动件的弹性变形； 4、 悬绳与悬绳器在垂直方向上的等效弹簧常数为定值；5、 油管全长锚定或不锚定，当油管不锚定时，考虑液体静负荷引起的油管柱弹性变形；6、 不考虑油管内液柱的振动；7、 假设抽油杆柱与油管同心；8、 油井是铅直的；9、 连杆与游梁之间的轴承座质量忽略；10、游梁与支架之间的轴承座质量忽略；11、所有构件为刚性构件；12、曲柄为等效构件；13、求等效力矩时不考虑全部

6、构件的质量；14、求等效转动惯量时考虑连杆质量；15、悬点载荷简化为F16、减速箱效率为0.866（1）抽油机运动规律建模根据抽油机型号，可以确定抽油机传动系统的机构尺寸、运动件质量、转动惯量、质心位置等参数。具体参数包括：R曲柄半径，m；P连杆长度，m；C游梁后臂长度，m；K基杆长度，m；A游梁前臂长度，m；I基杆的水平投影，m；曲柄平衡重滞后角，rad；JB1小皮带轮转动惯量，kg.m2；JB2大皮带轮转动惯量，kg.m2；JR1二级减速箱高速轴（输入轴）系转动惯量，kg.m2；JR2二级减速箱中间轴系转动惯量，kg.m2；JR3二级减速箱底速轴（输出轴）系转动惯量，kg.m2；iB皮带传

7、动比；iR1二级减速箱的第一级传动比；i R2二级减速箱的第二级传动比；mQ单曲柄质量，kg； rQ曲柄质心半径，m；MC单侧曲柄平衡块质量，kg； RC曲柄平衡块质心半径，m；JC单侧曲柄平衡块相对其质心的转动惯量，kg.m2；ML连杆质量，kg； JL连杆相对其质心的转动惯量，kg.m2；JB游梁组件相对其回转中心的转动惯量，kg.m2。 (3)油井参数油井参数包括油藏参数与生产参数。根据假设，和数学建模相关的几个参数包括：MR抽油杆柱质量，kg；抽油杆柱在液柱中的重量，N；油井动液面以上整个柱塞面积上的液柱重量，N；抽油杆柱在液体载荷作用下的静变形，m。各参数具体数值可取：R1.064m

8、；P3.026m；C2.82m；K=4.265m；A3.935m；I2.82m；0；JB12.0kg.m2；JB215.0kg.m2；JR13.0 kg.m2；JR28.0 kg.m2；JR315.0 kg.m2；iB=5.0；iR1=5.5；i R2=5.73；mQ=500kg； rQ=1.0m；MC=1000kg； RC=1.2m；JC=100kg.m2；ML=300kg； JL=300kg.m2；JB=1000kg.m2;MR=4000kg；=36000N；=30000N；=0.2m。游梁式抽油机示意图：（1） 各杆件角速度与角加速度图示K为初始位置 曲柄与初始位置夹角 连杆相对于初始位

9、置夹角摇杆相对于初始位置 将图中各杆件看成矢量，可以得到如下的矢量方程上述矢量方程可以用复变矢量表示为将上式两边对时间t求导，可得到各杆件的角速度为：将上式对时间求导，可得各杆件运动的角加速度： （2） 悬点位移、速度与角加速度悬点中程长度SPR: 悬点位移PR以下死点为位移零点，向上为位移的正方向，则任意时刻的悬点位移PR为： 另外以上死点为位移零点，向下为位移的正方向，则悬点位移为： 悬点加速度与加速度为：悬点向上运动为位移的正方向，则与为：扭矩因数（2）曲柄运动规律建模电机参数：电机型号额定功率/kw同步转速（r/m）满载转速（r/min）过载系数效率/%Y200l2-622100097

10、02.090.2（1）等效驱动力矩 因无实测的电动机机械特性曲线，所以根据电动机额定参数有如下计算：其中： （2）等效阻力矩段时（0，0.2）BC段时(0.2，SPR)CD段时(SPR,SPR-0.2)DA段时(SPR-0.2，0)(3)等效转动惯量其中各部件的动能：小皮带轮： 其中大皮带轮： 其中减速箱高速轴： 其中减速箱中间轴： 其中减速箱低速轴： 其中曲柄： 其中平衡块： 连杆：其中游梁组件： 其中悬点： 其中电动机转子： 其中（3） 曲柄扭矩AB段:M= (PR)-B- BC段：M= ()-B- CD段：M= -B-DA段：M= （-B）- = B= = (5)悬点载荷1,2曲线是抽油

11、机上冲程3,4曲线是抽油机下冲程抽油杆在液柱中的重量36000N抽油杆液面以上整个柱塞面积上的液柱重量30000N悬点冲程长度抽油杆在液体载荷作用下的静变形，0.2m坐标载荷：AB段 BC段CD段DA段 （6）系统运动微分方程其中：三 数学仿真方法（1）数学工具简介我们这次主要应用了MATLAB软件进行数学分析仿真。MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在

12、数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。在数值分析法上MATLAB具有巨大的优势，MATLAB的矩阵和向量操作功能是其他语言无法比拟的。在MATLAB环境下，数组的操作与数的操作一样简单，基本数据单元是不需要指定维数的，不需要说明数据类型的矩阵，而其数学表达式和运算规则与通常的习惯相同。MATLAB有着丰富的函数库，可以随时调用，大大方便了操作，在数值分析法上MATLAB具有巨大的优势。（2）龙格库塔法简介我们此次采用的数学方

13、法的是龙格库塔法，龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。龙格库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前而改进的欧拉公式将导数项取为

14、两端导数的平均。龙格-库塔方法的基本思想：在区间内多取几个点，将他们的斜率加权平均，作为导数的近似。龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。（4） 部分代码% 主程序solve1.m %求解微分方程clcclose allclear all% 最后修改于2016年6月12日% 4阶龙格库塔法微分方程求解% 声明全局变量global R P C K A I tau JB1 JB2 JR1 ., JR2 JR3 iB iR1 iR2 iMB yitaMB mQ rQ Mc Rc Jc ML JL JB MR PR0 PL0 k0 g ., kk n0

15、 nH PH yitaB w0 SH Sm MH;R=1.064;P=3.026;C=2.82;K=4.265;A=3.935;I=2.82;tau=0;%表示JB1=2;JB2=15;JR1=3;JR2=8;JR3=15;iB=5;iR1=5.5;iR2=5.73;iMB=iB*iR1*iR2;%电机到输出轴的传动比yitaMB=0.86;%电机到输出轴的效率mQ=500;rQ=1;%Mc=1000; %5500%Rc=1.2; %1.5Mc=1000; %调节Mc和Rc参数使Mef近似正弦Rc=1.04; Jc=100;ML=300;JL=300;JB=1000;MR=4000;PR0=3

16、6000;PL0=30000;k0 =0.2;%表示g = 9.8;%电机参数kk = 2;%表示kn0 = 1000;nH = 970;PH = 22;yitaB = 0.902;w0 = 2*pi*n0/60;SH = (n0-nH)/n0;Sm = SH*(kk+sqrt(kk2-1);MH = 9550*PH/nH;% 外层求得Je的导函数global diffJe_thetai;global diffJe THETA;syms THETA;a = asin(I/K);theta0 = pi/2-asin(P2+K2-C2)./(2*P*K)-a;THETA1=theta0+THETA

17、;THETA2=2*pi-THETA1+a;Ll = sqrt(R2+K2-2*R*K*cos(THETA2);BETA = asin(R*sin(THETA2)/Ll);THETA3 = acos(P2+Ll2-C2)/(2*P*Ll) - BETA;THETA4 = acos(P2-Ll2-C2)/(2*P*Ll) - BETA;%关于theta的自变量DJe(THETA) = JB1*iMB2+(JB2+JR1)*(iR1*iR2)2+JR2*iR22+JR3+2*mQ*rQ2+ ., 2*mQ*rQ2/3+2*Jc+2*Mc*Rc2+JL*(R*sin(THETA4-THETA2)/(

18、P*sin(THETA3-THETA4)2+ ., ML*( JL*R2*(1+ (sin(THETA4-THETA2)2/(sin(THETA3-THETA4)2/4 + sin(THETA4-THETA2)*cos(THETA2-THETA3)/sin(THETA3-THETA4) ) )+ ., JB*(R*sin(THETA4-THETA2)/P/sin(THETA3-THETA4)2+MR*(A*R*sin(THETA3-THETA2)/P/sin(THETA3-THETA4)2;diffJe = diff(DJe(THETA);h = 0.01; %步长t = 0:h:20; %时

19、域RK_theta(1)=0; %初始转角RK_w(1)=0; %初始角速度for i = 1:length(t)-1 K1 = RK_w(i); diffJe_thetai = subs(diffJe,RK_theta(i); L1=rightf_sys1(t(i),RK_theta(i),RK_w(i); K2=RK_w(i)+0.5*h*L1; diffJe_thetai = subs(diffJe,RK_theta(i)+0.5*h*K1); L2=rightf_sys1(t(i)+0.5*h,RK_theta(i)+0.5*h*K1,RK_w(i)+0.5*h*L1); K3=RK_

20、w(i)+0.5*h*L2; diffJe_thetai = subs(diffJe,RK_theta(i)+0.5*h*K2); L3=rightf_sys1(t(i)+0.5*h,RK_theta(i)+0.5*h*K2,RK_w(i)+0.5*h*L2); K4=RK_w(i)+h*L3; diffJe_thetai = subs(diffJe,RK_theta(i)+h*K3); L4=rightf_sys1(t(i)+h,RK_theta(i)+h*K3,RK_w(i)+h*L3); %RK_theta(i+1)=mod(RK_theta(i)+1/6*h*(K1+2*K2+2*K3

21、+K4),(2*pi); RK_theta(i+1)=RK_theta(i)+1/6*h*(K1+2*K2+2*K3+K4); RK_w(i+1) =RK_w(i)+1/6*h*(L1+2*L2+2*L3+L4);end for i = 1:length(t) RK_theta(i)=mod(RK_theta(i),(2*pi); end%保存数值解到datasave data RK_theta RK_w tplot(t,RK_theta,r);grid onxlabel(时间t/s);ylabel(曲柄转角/rad);legend(t);figureplot(t(1:150),RK_w(1:

22、150),b);xlabel(时间t/s);ylabel(曲柄角速度/(rad/s);grid onlegend(t);四 系统动力特性仿真分析（1） 电机相关曲线（2） 悬点图（3） 杆件运动规律图（4）等效阻力矩(5)曲柄扭矩曲线图（6）等效转动惯量曲线（7）静力示功图五 收获与感想通过这次课程设计，我们深深的意识到自己还有很多知识需要学习，自己也还有很多地方有待提高。从之前对系统动力学的初步了解，到项目讨论完成后，我们对动力学有了进一步的认识和了解，不仅仅使我们复习了课程的内容，也加深了对线性微分方程组解法的了解，对以后解方程有了更深入的了解。在项目讨论的过程中，我们查阅了很多关于算法的

23、书籍及文献，也锻炼了我们查阅资料的能力。在这次项目过程中，我们用到了matlab软件，刚接触这个软件，感觉很陌生，无从下手，经过这些天的学习与应用，对其有了进一步的了解，这巩固了我们对知识的了解，也在一定程度上锻炼了自己对matlab软件的开发能力。从团队的角度来说，大家都为这个项目付出了很多，大家在一起互相讨论，激智与集智于一体，互相学习彼此的优点。六 参考文献1.抽油机设计计算与计算机实现 董世民 张士军 编著2.机械系统动力学 董世民 编著3.机械原理 安子军 编著4. 理论力学 杜国君 编著七 附录（程序代码）% 主程序solve1.m %求解微分方程clcclose allclear

24、 all% 最后修改于2016年6月12日% 4阶龙格库塔法微分方程求解% 声明全局变量global R P C K A I tau JB1 JB2 JR1 ., JR2 JR3 iB iR1 iR2 iMB yitaMB mQ rQ Mc Rc Jc ML JL JB MR PR0 PL0 k0 g ., kk n0 nH PH yitaB w0 SH Sm MH;R=1.064;P=3.026;C=2.82;K=4.265;A=3.935;I=2.82;tau=0;%表示JB1=2;JB2=15;JR1=3;JR2=8;JR3=15;iB=5;iR1=5.5;iR2=5.73;iMB=iB

25、*iR1*iR2;%电机到输出轴的传动比yitaMB=0.86;%电机到输出轴的效率mQ=500;rQ=1;%Mc=1000; %5500%Rc=1.2; %1.5Mc=1000; %调节Mc和Rc参数使Mef近似正弦Rc=1.04; Jc=100;ML=300;JL=300;JB=1000;MR=4000;PR0=36000;PL0=30000;k0 =0.2;%表示g = 9.8;%电机参数kk = 2;%表示kn0 = 1000;nH = 970;PH = 22;yitaB = 0.902;w0 = 2*pi*n0/60;SH = (n0-nH)/n0;Sm = SH*(kk+sqrt(

26、kk2-1);MH = 9550*PH/nH;% 外层求得Je的导函数global diffJe_thetai;global diffJe THETA;syms THETA;a = asin(I/K);theta0 = pi/2-asin(P2+K2-C2)./(2*P*K)-a;THETA1=theta0+THETA;THETA2=2*pi-THETA1+a;Ll = sqrt(R2+K2-2*R*K*cos(THETA2);BETA = asin(R*sin(THETA2)/Ll);THETA3 = acos(P2+Ll2-C2)/(2*P*Ll) - BETA;THETA4 = acos

27、(P2-Ll2-C2)/(2*P*Ll) - BETA;%关于theta的自变量DJe(THETA) = JB1*iMB2+(JB2+JR1)*(iR1*iR2)2+JR2*iR22+JR3+2*mQ*rQ2+ ., 2*mQ*rQ2/3+2*Jc+2*Mc*Rc2+JL*(R*sin(THETA4-THETA2)/(P*sin(THETA3-THETA4)2+ ., ML*( JL*R2*(1+ (sin(THETA4-THETA2)2/(sin(THETA3-THETA4)2/4 + sin(THETA4-THETA2)*cos(THETA2-THETA3)/sin(THETA3-THET

28、A4) ) )+ ., JB*(R*sin(THETA4-THETA2)/P/sin(THETA3-THETA4)2+MR*(A*R*sin(THETA3-THETA2)/P/sin(THETA3-THETA4)2;diffJe = diff(DJe(THETA);h = 0.01; %步长t = 0:h:20; %时域RK_theta(1)=0; %初始转角RK_w(1)=0; %初始角速度for i = 1:length(t)-1 K1 = RK_w(i); diffJe_thetai = subs(diffJe,RK_theta(i); L1=rightf_sys1(t(i),RK_th

29、eta(i),RK_w(i); K2=RK_w(i)+0.5*h*L1; diffJe_thetai = subs(diffJe,RK_theta(i)+0.5*h*K1); L2=rightf_sys1(t(i)+0.5*h,RK_theta(i)+0.5*h*K1,RK_w(i)+0.5*h*L1); K3=RK_w(i)+0.5*h*L2; diffJe_thetai = subs(diffJe,RK_theta(i)+0.5*h*K2); L3=rightf_sys1(t(i)+0.5*h,RK_theta(i)+0.5*h*K2,RK_w(i)+0.5*h*L2); K4=RK_w(

30、i)+h*L3; diffJe_thetai = subs(diffJe,RK_theta(i)+h*K3); L4=rightf_sys1(t(i)+h,RK_theta(i)+h*K3,RK_w(i)+h*L3); %RK_theta(i+1)=mod(RK_theta(i)+1/6*h*(K1+2*K2+2*K3+K4),(2*pi); RK_theta(i+1)=RK_theta(i)+1/6*h*(K1+2*K2+2*K3+K4); RK_w(i+1) =RK_w(i)+1/6*h*(L1+2*L2+2*L3+L4);end for i = 1:length(t) RK_theta(

31、i)=mod(RK_theta(i),(2*pi); end%保存数值解到datasave data RK_theta RK_w tplot(t,RK_theta,r);grid onxlabel(时间t/s);ylabel(曲柄转角/rad);legend(t);figureplot(t(1:150),RK_w(1:150),b);xlabel(时间t/s);ylabel(曲柄角速度/(rad/s);grid onlegend(t);子函数区rightf_sys1.mfunction alpha = rightf_sys1(t,theta,w)% 参数定义 使用全局变量global R P C K A I tau JB1 JB2 JR1 ., JR2 JR3 iB iR1 iR2 iMB yitaMB mQ rQ Mc Rc Jc ML JL JB MR PR0 PL0