valve_motor_pinion.m

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

%valve_moter_pinion.m
%differential equation
function dxdt=valve_motor_pinion(t,x)
%dxdt=zeros(8,1);
%x=zeros(8,1);
%ui=zeros(4,1);
ui=[10;-10;10;-10];
%------------------------
i1=44.8;
i2=4.21;
%--------------------油路1的参数
Jm1=200;        %马达1的转动惯量
Ip1=4000;       %小齿轮1的转动惯量
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=600000;    %马达1的负载力矩
Tmf1=10000;   %马达1的摩擦力矩
Tf1=100000;   %齿轮副1的摩擦力矩
Te1=(i1*i2-1)*TL1-Tmf1-Tf1*i2;    %等效力矩
Dm1=159.23;   %马达排量
beta1=1.4e+9;   %油体积弹性模量
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=100000;     %马达2的负载力矩
Tmf2=10000;     %马达2的摩擦力矩
Tf2=100000;      %齿轮副2的摩擦力矩
Te2=(i1*i2-1)*TL2-Tmf2-Tf2*i2;    %等效力矩
Dm2=159.23;      %马达2排量
beta2=1.4e+9;      %马达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的回油口压力
%----------------------------------
%differential equation for valve,motor,pinion
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));