|
%初始条件
theta1=linspace(0,360,100);%单位度
theta1=theta1*pi/180;%转换为弧度制
W1=52*pi/30;%角速度 单位rad/s
H=0.33;%行程 单位m
Lo2a=0.099;%O2A的长度 单位m
Lo4b=0.599;%O4B的长度 单位m
Lbc=0.198;%Bc的长度 单位m
Lo2o4=0.360;%O2O4的长度 单位m
Lo4d=0.576;
Z=pi/180;%角度与弧度之间的转换
%求解Lo4a、Theta4、Theta5和ldc四个变量
Lo4a=((Lo2o4)^2+(Lo2a)^2-2*Lo2o4*Lo2a*cos(theta1+pi/2)).^0.5;%求出O4A的值
for i=1:100%求解角度theta4、Theta5和Ldc的长度
theta4(i)=acos(Lo2a*cos(theta1(i))/Lo4a(i));
theta5(i)=asin((Lo4d-Lo4b*sin(theta4(i)))/Lbc); Ldc(i)=Lo4b*cos(theta4(i))+Lbc*cos(theta5(i));
end%求解完成
%求解完成
%求解VLo4a、W4、W5和Vc四个变量
for i=1:100
J=inv([cos(theta4(i)),-Lo4a(i)*sin(theta4(i)),0,0; sin(theta4(i)),Lo4a(i)*cos(theta4(i)),0,0; 0,-Lo4b*sin(theta4(i)),-Lbc*sin(theta5(i)),-1; 0,Lo4b*cos(theta4(i)),Lbc*cos(theta5(i)),0]); K=J*W1*[-Lo2a*sin(theta1(i));Lo2a*cos(theta1(i));0;0];
VLo4a(i)=K(1);
W4(i)=K(2);
W5(i)=K(3);
Vc(i)=K(4);
end%求解完成
%求解aLo4a、a4、a5、ac四个变量
for i=1:100
J=inv([cos(theta4(i)),-Lo4a(i)*sin(theta4(i)),0,0; sin(theta4(i)),Lo4a(i)*cos(theta4(i)),0,0; 0,-Lo4b*sin(theta4(i)),-Lbc*sin(theta5(i)),-1; 0,Lo4b*cos(theta4(i)),Lbc*cos(theta5(i)),0]);
P=W1*W1*[-Lo2a*cos(theta1(i));-Lo2a*sin(theta1(i));0;0];
M=[-W4(i)*sin(theta4(i)),-VLo4a(i)*sin(theta4(i))-Lo4a(i)*W4(i)*cos(theta4(i)),0,0;
W4(i)*cos(theta4(i)),VLo4a(i)*cos(theta4(i))-Lo4a(i)*W4(i)*sin(theta4(i)),0,0;
0,-Lo4b*W4(i)*cos(theta4(i)),-Lbc*W5(i)*cos(theta5(i)),0; 0,-Lo4b*W4(i)*sin(theta4(i)),-Lbc*W5(i)*sin(theta5(i)),0]; N=[VLo4a(i);W4(i);W5(i);Vc(i)];
K=J*(-M*N+P);
aLo4a(i)=K(1);
a4(i)=K(2);
a5(i)=K(3);
ac(i)=K(4);
end%求解完成
%动态静力分析
%初始条件
G4=220; G6=800; Js4=1.2; P=9000; Lo4s4=0.5*Lo4b;Xs6=200;Ys6=50;Yp=80;Yo4=0;Xo4=0;Yo2=0.36;Xo2=0;
%给切削阻力赋值
for i=1:100
if((abs(Ldc(1)-Ldc(i))>0.05*H&&abs(Ldc(1)-Ldc(i))<0.95*H)&&(theta1(i)<pi))
P(i)=9000; else
P(i)=0; end
end%赋值完成
F16(i)=(G6/10)*ac(i);
Fr6(i)=(F16(i)*Ys6+P(i)*Yp+G6*Xs6)/Xs6;
Ya(i)=Lo4a(i)*sin(theta4(i));
Xa(i)=Lo4a(i)*sin(theta4(i));
Ys4(i)=Lo4s4*sin(theta4(i));
Xs4(i)=Lo4s4*cos(theta4(i));
Yb(i)=Lo4b*sin(theta4(i));
Xb(i)=Lo4b*cos(theta4(i));
F4x(i)=-(G4/10)*(a4(i)*Lo4s4*sin(theta4(i)));
F4y(i)=-(G4/10)*(a4(i)*Lo4s4*cos(theta4(i)));
%求解平衡力矩
J4=Js4+(G4/10)*(0.5*Lo4b)*(0.5*Lo4b);%导杆对点O4的转动惯量
M4(i)=-J4*a4(i);
J=inv([1,0;
0,1]);
K=J*[ F16(i)+P(i);-Fr6(i)+G6]
F56x(i)=K(1);F56y(i)=K(2);
F45x(i)= F56x(i);
F45y(i)= F56y(i);
Y1= ' F43x(i)=F23x(i)) ' ;
Y2= 'F43y(i)=F23y(i) ' ;
Y3= 'F43x(i)*abs(Lo2a*sin(theta1(i)))+F43y(i)*abs(Lo2a*cos(theta1(i)))=abs(Lo2a*sin(theta1(i)))*F23x(i)+ abs(Lo2a*cos(theta1(i)))* F23y(i) ' ;
Y4= ' F14x(i)+F4x(i)=F45X(i)+F43x(i)' ;
Y5= 'F4y(i)+ F14y(i)=(-F4y(i)F43y(i))+F45y(i) ' ;
Y6= ' abs(Ya(i))* F43x(i)+abs(Xa(i))*F34y+abs(Yb(i))*F54x+abs(Xb(i))*F54y+M4+G4*abs(Xs4)=0 ' ;
Y7= ' F12x(i)=F23x(i) ' ;
Y8= 'F12y(i)=F23y(i) ' ;
[K1,K2,K3,K4,K5,K6,K7,K8]=solve(Y1,Y2,Y3,Y4,Y4,Y5,Y6,Y7,Y8,' F23x(i), F23y(i), F43x(i), F43y(i), F14x(i), F14y(i),F12x(i), F12y(i) ');
F23x(i)=K1;F23y(i)=K2;F43x(i)=K3;F43y(i)=K4;
F14x(i)=K5;F14y(i)=K6;F12x(i)=K7;F12y(i)=K8;
F56(i)=sqrt(F56x(i)*F56x(i)+F56y(i)*F56y(i));
F45(i)=sqrt(F45x(i)*F45x(i)+F45x(i)*F45y(i));
F23(i)=sqrt(F23x(i)*F23x(i)+F23y(i)*F23y(i));
F43(i)=sqrt(F43x(i)*F43x(i)+F43y(i)*F43y(i));
F14(i)=sqrt(F14x(i)*F14x(i)+F14y(i)*F14y(i));
F12(i)=sqrt(F12x(i)*F12x(i)+F12y(i)*F12y(i));
u(i)=F23x(i)/F23y(i);
t(i)=atan(u(i));
M(i)=F32(i)*Lo2a*sin(pi*3/2-theta1(i)-t(i));
for i=1:100
Ekk(i)=((G6/10)*Vc(i)*Vc(i))/2;%计算总动能
end
dEkk(1)=Ekk(1)-Ekk(100);%动能的改变量
for i=2:100
dEkk(i)=Ekk(i)-Ekk(i-1);%动能的改变量
end
for i=1:100
MM(i)=(dEkk(i)+P(i)*abs(Vc(i)))/W1;%求平衡力矩
End
for i= 1:100
theta1(i)=(theta1(i)*180)/pi;
theta4(i)=(theta4(i)*180)/pi;
theta5(i)=(theta5(i)*180)/pi;
%画图
%画运动图
figure(1);
plot(theta1,theta4,'r',theta1,theta5,'b');
grid on;
xlabel('角度（度）');
ylabel('theta4、theta5(度)');
title('角度Theta4、theta5');
axis([0,360,0,360]);
figure(2);
plot(theta1,Ldc);grid on;
xlabel('角度（度）');
ylabel('Ldc(m)');
title('刨头6位移');
axis('auto');
figure(3);
plot(theta1,W4,'r',theta1,W5,'b');grid on;
xlabel('角度（度）');
ylabel('W4、W5(rad/s)');
title('角速度W4、W5');
axis([0 , 360,-5,3]);
figure(4);
plot(theta1,Vc);grid on;
xlabel('角度（度）');
ylabel('Vc(m/s)');
title('刨头6的速度');
axis([0 , 360,-5,3]);
figure(5);
plot(theta1,a4,'r',theta1,a5,'b');grid on;
xlabel('角度（度）');
ylabel('a4、a5(rad/s/s)');
title('角度加速度a4、a5');
axis([0,360,-80,80]);
figure(6);
plot(theta1,ac,'r');grid on;
xlabel('theta1角度（度）');
ylabel('ac(m/s/s)');
title('刨头加速度ac');
axis([0,360,-80,80]);
%运动图画完
%画反力图
figure(7);
plot(theta1,F56,'r',theta1,F45,'b',theta1,F23,'g');grid on;
xlabel(' theta1角度（度）');
ylabel('F56、F45、F23(rad/s/s)');
title('运动副反力F56、F45、F23');
axis('auto');
figure(8);
plot(theta1,F43,'r',theta1,F14,'b',theta1,F12,'g');grid on;
xlabel(' theta1角度（度）');
ylabel('F43、F14、F12(rad/s/s)');
title('运动副反力F43、F14、F12');
axis('auto');
figure(9);
plotyy(theta1,P,theta1,Ldc);
xlabel(' theta1角度（度）');
ylabel('P');
axis([theta1(1) ,theta1(100),-50,1400]);
title('切削阻力P与位移Ldc');grid on;
figure(10);
plot(theta1,MM,'r');grid on;
xlabel(' theta1角度（度）');
ylabel('MM(N·m)');
figure(11);
plot(theta1,MM,'r');grid on;
xlabel(' theta1角度（度）');
ylabel('M(N·m)');
title(' 平衡力矩')
上述是所有程序
J=inv([1,0;
0,1]);
K=J*[ F16(i)+P(i);-Fr6(i)+G6]
F56x(i)=K(1);F56y(i)=K(2);
F45x(i)= F56x(i);
F45y(i)= F56y(i);
Y1= ' F43x(i)=F23x(i)) ' ;
Y2= 'F43y(i)=F23y(i) ' ;
Y3= 'F43x(i)* abs(Lo2a*sin(theta1(i)))+F43y(i)*abs(Lo2a*cos(theta1(i)))=abs(Lo2a*sin(theta1(i)))* F23x(i)+ abs(Lo2a*cos(theta1(i)))* F23y(i) ' ;
Y4= ' F14x(i)+F4x(i)=F45X(i)+F43x(i) ' ;
Y5= 'F4y(i)+ F14y(i)= (-F4y(i)F43y(i))+F45y(i) ' ;
Y6= ' abs(Ya(i))* F43x(i)+ abs(Xa(i))*F34y+abs(Yb(i))*F54x+abs(Xb(i))*F54y+M4+G4*abs(Xs4)=0 ' ;Y7= ' F12x(i)=F23x(i) ' ;
Y8= 'F12y(i)=F23y(i) ' ;
[K1,K2,K3,K4,K5,K6,K7,K8]=solve(Y1,Y2,Y3,Y4,Y4,Y5,Y6,Y7,Y8,' F23x(i), F23y(i), F43x(i), F43y(i), F14x(i), F14y(i), F12x(i), F12y(i) ');
主要是这段程序的错误，矩阵为奇异工作精度，解方程组时出现一堆错误 请大神们帮我修改下程序
| 
| Archiver|手机版|小黑屋| 碧波制图网
Published by Stonespider 
