end
for i=(N/2+1):N;
DJ1DF1(i)=(RJ1(i+1)-RJ1(i))./(XTA(i+1)-XTA(i))./2.0+(RJ1(i)-RJ1(i-1))./(XTA(i)-XTA(i-1))./2.0;
end
a=RJ1;
b=DJ1DF1;
c=F1+FM1;
dy(1)=y(2);
dy(2)=-b*(y(2))^2/2/a+c/a;
全满弹时弹仓的运动微分方程:
%弹筒的运动微分方程
function dy=dantong(t,y)
dy=zeros(2,1);
%求各时刻,各转角下的等效转动惯量及其对转角的导数
%求各时刻,个转角下的等效力矩(所有弹筒的重力产生的等效力矩,所有弹筒摩擦力产生的等效力矩)
%RJ-不同时刻各弹筒的等效转动惯量
%EJ-不同时刻所有弹筒的等效转动惯量
%F-不同时刻各(每个)弹筒重力产生的等效力矩
%EF-不同时刻所有弹筒的重力产生的等效力矩
%FM-不同时刻各(每个)弹筒摩擦力产生的等效力矩
%EFMDT-不同时刻所有弹筒摩擦力产生的等效力矩
OMGA=100.0/180.0*pi;
DT=0.5e-3;
N=pi/2.0/OMGA/DT;
for IT=1:(N+1)
EJ(IT)=0.0;
EF(IT)=0.0;
EFMDT(IT)=0.0;
T=(IT-1)*DT;
XTA(IT)=T*OMGA;
RM1=68.8398; %RM1--弹筒1的质量,由计算机测量获得
RJ1=0.294227; %RJ1--弹筒1的转动惯量,由计算机测量获得
U=0.15; %U--摩擦系数
OB=131e-3; %参数(1)
BA=sqrt(2*131e-3^2); %参数(2)
OB1=131e-3; %参数(3)
CA=sqrt(72e-3^2+92.63e-3^2); %参数(4)
CAB=atan(72e-3/92.63e-3); %参数(5)
R=pi/180.0;
BB1=sqrt(OB^2+OB^2-2.0*OB^2*cos(XTA)); %(1-1)
BB1O=(pi-XTA)/2.0; %(1-2)
BB1A=pi/2+BB1O; %(1-3)
B1AB=asin(BB1./BA.*sin(BB1A)); %(1-4)
ABB1=pi-BB1A-B1AB; %(1-5)
B1A=sqrt(BB1.^2+BA.^2-2.*BB1.*BA.*cos(ABB1)); %(1-6)
OA=sqrt(B1A.^2+OB1.^2); %(1-7)
S1=BA-B1A; %(1-8)
OBA=acos((BA.^2+OB.^2-OA.^2)./(2.0.*BA.*OB)); %(1-9)
ALF=OBA; %(1-10)
VB=OB.*OMGA; %(1-11)
VZ=VB.*sin(XTA)./sin(pi-XTA-ALF); %(1-12)
OMGAZ=VZ./BA; %(1-13)
VP1=VB.*sin(ALF)./sin(pi-XTA-ALF); %(1-14)
CAB1=CAB-B1AB; %(1-15)
BTA=CAB1; %(1-16)
GMA=pi-(pi./2.0-BTA); %(1-17)
VC=sqrt(VP1.^2+(CA.*OMGAZ).^2-2.0.*VP1.*CA.*OMGAZ.*cos(GMA)); %(1-18)
ALFG=acos(CA.*OMGAZ.*cos(BTA)./VC); %(1-19)
RJ1=RM1.*(VC./OMGA).^2+RJ1.*(OMGAZ./OMGA).^2;
F1=RM1.*9.8.*VC.*cos(ALFG)./OMGA;
FM1=RM1.*9.8/2.0.*U.*VP1/OMGA.*cos(pi);
RM2=68.8398; %RM2--弹筒2的质量,由计算机测量获得
RJ2=0.294227; %RJ2--弹筒2的转动惯量,由计算机测量获得
U=1.5; %U--摩擦系数
OC=92.63e-3+72.0e-3;
C1OB1=pi/4.0;
OMGAZ=OMGA;
VC=OC.*OMGA;
ALFG=abs(C1OB1-XTA);
RJ2=RM2.*(VC./OMGA).^2+RJ2.*(OMGAZ./OMGA).^2;
F2=RM2.*9.8.*VC.*cos(ALFG)./OMGA;
FM2=0;
%计算主动组合链轮匀速转动时弹筒3的运动规律:认为弹筒3做平面运动,求弹筒3的平均速度、转动速度及其质心的速度。
RM3=68.8398;
RJ3=0.294227;
U=0.15;
AO=131e-3;
B1O=131e-3;
AB=sqrt(2.0*131.0e-3^2);
ABC=atan(72.0e-3/92.63e-3);
CB=sqrt(72.0D-3^2+92.63D-3^2);
AB1=sqrt(AO.^2+B1O.^2-2.*AO.*B1O.*cos(pi/2-XTA));
B1AO=(pi-(pi/2-XTA))/2;
OB1A=(pi-(pi/2-XTA))/2;
BB1A=pi/2+OB1A;
ABB1=asin(AB1./AB.*sin(BB1A));
BAB1=pi-BB1A-ABB1;
B1B=sqrt(AB.^2+AB1.^2-2.0.*AB.*AB1.*cos(BAB1));
S3=B1B;
OB=sqrt(B1O.^2+B1B.^2);
BAO=acos((AO.^2+AB.^2-OB.^2)./(2.0.*AO.*AB));