valve_motor_pinion.asv

自己编的四节龙格琨塔程序

%valve_moter_pinion.m
%differention equation
function dxdt=valve_motor_pinion(t,x,ui)
i1=44.8;
i2=4.21;
%--------------------油路1的参数
Jm1=200;%马达的转动惯量
Ip1=4000;%小齿轮的转动惯量
Je1=Jm1*i1+Ip1*i2;%等效惯量
Bm1=1e-3;%马达1的阻尼
Rp1=0.300;%小齿轮1的分度圆半径
Cm1=0.005;%马达泄漏系数
Cn1=4000;%齿轮副1的啮合阻尼
Rg=0.600;%大齿轮的分度圆半径
Be1=Bm1*i1+Rp1*Cn1*(Rp1*i2-Rg);%等效阻尼
G1=5e+6;%马达1的输出轴刚度
Kn1=9e+4;%齿轮副1的啮合刚度
Ge1=G1*i1+Rp1*Kn1*(Rp1*i2-Rg);%等效刚度
TL1=60000;%马达1的负载力矩
Tmf1=10000;%马达1的摩擦力矩
Tf1=100000;%齿轮副1的摩擦力矩
Te1=(i1*i2-1)*TL1-Tmf1-Tf1*i2;%等效力矩
Dm1=159.23;%马达排量
beta1=2e+5;%油体积弹性模量
V01=500;%马达单腔容积
Cd=0.62;%先导阀阀芯的流量系数
omega=1.63e-2;%先导阀阀口的面积梯度
row=874;%油液密度
Kap1=1333;%感应式位移传感器1的增益
Kap2=1333;%感应式位移传感器2的增益
Pp1=3.15e+7;%阀1的输入油压力
Pp2=3.15e+7;%阀2的输入油压力
PT1=6e+5;%阀1的回油口压力
PT2=6e+5;%阀2的回油口压力
%------------------------------油路2的参数
Jm2=200;%马达2的转动惯量
Ip2=4000;%小齿轮2的转动惯量
Je2=Jm2*i1+Ip2*i2;%等效惯量
Bm2=1e-3;%马达2的阻尼
Rp2=0.300;%小齿轮2的分度圆半径
Cm2=0.005;%马达2泄漏系数
Cn2=4000;%齿轮副2的啮合阻尼
Rg=0.6;%大齿轮的分度圆半径
Be2=Bm2*i1+Rp2*Cn2*(Rp2*i2-Rg);%等效阻尼
G2=5e+6;%马达2的输出轴刚度
Kn2=9e+4;%齿轮副2的啮合刚度
Ge2=G2*i1+Rp2*Kn2*(Rp2*i2-Rg);%等效刚度
TL2=60000;%马达2的负载力矩
Tmf2=10000;%马达2的摩擦力矩
Tf2=100000;%齿轮副2的摩擦力矩
Te2=(i1*i2-1)*TL2-Tmf2-Tf2*i2;%等效力矩
Dm2=159.23;%马达2排量
beta2=2e+5;%马达2油体积弹性模量
V02=500;%马达2单腔容积
Cd=0.62;%先导阀阀芯的流量系数
omega=1.63e-2;%先导阀阀口的面积梯度
row=874;%油液密度
Kap3=1333;%感应式位移传感器3的增益
Kap4=1333;%感应式位移传感器4的增益
Pp3=3.15e+7;%阀3的输入油压力
Pp4=3.15e+7;%阀4的输入油压力
PT3=6e+5;%阀3的回油口压力
PT4=6e+5;%阀4的回油口压力
%----------------------------------
%difference equation for valve,motor,pinion
dxdt=zeros(8,1);
x=zeros(8,1);
ui=zeros(4,1);
dxdt(1)=x(2);
dxdt(2)=(1/Je1)*(-Ge1*x(1)-Be1*x(2)+Dm1*x(3)-Dm1*x(4)+Te1);
dxdt(3)=(beta1/V01)*(-Dm1*i1*x(2)-Cm1*x(3)+Cm1*x(4)+(Cd*omega*ui(1)/Kap1)*sqrt(2/row)*sqrt((Pp1-x(3))/2+sign(ui(1))*(Pp1-x(3))/2+(x(3)-PT1)/2-sign(ui(1))*(x(3)-PT1)/2));
dxdt(4)=(beta1/V01)*(Dm1*i1*x(2)+Cm1*x(3)-Cm1*x(4)-(Cd*omega*ui(2)/Kap2)*sqrt(2/row)*sqrt((Pp2-x(4))/2+sign(ui(2))*(Pp2-x(4))/2+(x(4)-PT2)/2-sign(ui(2))*(x(4)-PT2)/2));
dxdt(5)=x(6);
dxdt(6)=(1/Je2)*(-Ge2*x(5)-Be2*x(6)+Dm2*x(7)-Dm2*x(8)+Te2);
dxdt(7)=(beta2/V02)*(-Dm2*i1*x(6)-Cm2*x(7)+Cm2*x(8)+(Cd*omega*ui(3)/Kap3)*sqrt(2/row)*sqrt((Pp3-x(7))/2+sign(ui(3))*(Pp3-x(7))/2+(x(7)-PT3)/2-sign(ui(3))*(x(7)-PT3)/2));
dxdt(8)=(beta2/V02)*(Dm2*i1*x(6)+Cm2*x(7)-Cm2*x(8)-(Cd*omega*ui(4)/Kap3)*sqrt(2/row)*sqrt((Pp4-x(8))/2+sign(ui(4))*(Pp4-x(8))/2+(x(8)-PT4)/2-sign(ui(4))*(x(8)-PT4)/2));