codice:
Sfunction mocen
% ****************************************************
% programma per la simulazione della cinematica del
% manovellismo ordinario centrato
% ****************************************************
% Inserimento dati
stato = 1;
while stato == 1
a = input('Assegna raggio manovella in [m]: ');
b = input('Assegna lunghezza biella in [m]: ');
if a < b
stato = 0;
else
disp('La manovella deve essere più corta della biella')
end
end
alfa = [1:0.1:360].*(pi/180);
alfap = input('Assegna velocita'' angolare manovella [rad/s]: ');
alfapp = 0;
% Analisi di posizione
disp('Risoluzione analisi di posizione')
c = a*cos(alfa) + b*sqrt(1-((a*sin(alfa))./b).^2);
beta = asin(-(a*sin(alfa))./b);
%% PRIMA E SECONDA APPROSSIMAZIONE
c1_approx=a*cos(alfa) + b; %Buona approx solo se lambda è molto minore di 1
c2_approx=a*cos(alfa) + b + a*a/b/4*(cos(2*alfa)-1) ;
%% PIEDE DI BIELLA
xb=a+b-c;
xb1_approx=a+b-c1_approx;
xb2_approx=a+b-c2_approx;
% Analisi di velocità
disp('Risoluzione analisi di velocità')
for i = 1:length(alfa)
A = [1 b*sin(beta(i));
0 -b*cos(beta(i))];
B = [-alfap*a*sin(alfa(i));
alfap*a*cos(alfa(i))];
xp = inv(A)*B;
cp(i) = xp(1);
betap(i) = xp(2);
end
% cp = -alfap*a*(sin(alfa) - cos(alfa).*tan(beta));
% betap = -alfap*((a*cos(alfa))./(b*cos(beta)));
%% PRIMA E SECONDA APPROSSIMAZIONE
cp1_approx=a*alfap * ( -sin(alfa) );
cp2_approx=a*alfap * ( -sin(alfa) -a/b/2*sin(2*alfa) );
% Analisi di accelerazione
disp('Risoluzione analisi di accelerazione')
for i = 1:length(alfa)
A = [1 b*sin(beta(i));
0 -b*cos(beta(i))];
B = -[ alfapp*a*sin(alfa(i)) + alfap^2*a*cos(alfa(i)) + betap(i)^2*b*cos(beta(i));
-alfapp*a*cos(alfa(i)) + alfap^2*a*sin(alfa(i)) + betap(i)^2*b*sin(beta(i))];
xpp = inv(A)*B;
cpp(i) = xpp(1);
betapp(i) = xpp(2);
end
%% PRIMA E SECONDA APPROSSIMAZIONE
cpp1_approx=a* alfap^2 * ( -cos(alfa) );
cpp2_approx=a* alfap^2 * ( -cos(alfa) -a/b*cos(2*alfa) );
% Visualizzazione risultati
stato = 1;
while stato == 1
k = menu('Visualizzazione risultati', ...
'Grafici posizione', ...
'Grafici velocita''', ...
'Grafici accelerazione', ...
'Animazione del moto', ...
'C.I.R.',...
'Fine');
if k == 1
figure(1)
subplot(211)
%plot(alfa,c)
plot(alfa,c,alfa,c1_approx,alfa,c2_approx)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Posizione piede di biella [m]'))
legend('Esatta','Prima approx','Seconda approx')
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
subplot(212)
plot(alfa,beta)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Rotazione biella [rad]\n'))
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
elseif k == 2
figure(2)
subplot(211)
%plot(alfa,cp)
plot(alfa,cp,alfa,cp1_approx,alfa,cp2_approx)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Velocita'' piede di biella [m/s]'))
legend('Esatta','Prima approx','Seconda approx')
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
subplot(212)
plot(alfa,betap)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Velocita'' angolare biella [rad/s]\n'))
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
figure(20)
plot(a-xb,-cp./alfap,(a-xb1_approx),-cp1_approx./alfap,(a-xb2_approx),-cp2_approx./alfap)
xlabel('r-x_B')
ylabel('v_B/\omega')
legend('Esatta','Prima approx','Seconda approx')
grid
axis equal
text(a*1.1,0,'PMS')
text(-a*1.25,0,'PMI')
elseif k == 3
figure(3)
subplot(211)
%plot(alfa,cpp)
plot(alfa,cpp,alfa,cpp1_approx,alfa,cpp2_approx)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Accelerazione piede di biella [m/s^2]'))
legend('Esatta','Prima approx','Seconda approx')
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
subplot(212)
plot(alfa,betapp)
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('Accelerazione angolare biella [rad/s^2]\n'))
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
figure(20)
plot((a-xb),-cpp./alfap^2,(a-xb1_approx),-cpp1_approx./alfap^2,(a-xb2_approx),-cpp2_approx./alfap^2)
xlabel('r-x')
ylabel('a_B/\omega^2')
legend('Esatta','Prima approx','Seconda approx')
grid
axis equal
text(a*1.1,0,'PMS')
text(-a*1.25,0,'PMI')
elseif k == 4
figure(4)
set(gcf,'DoubleBuffer','on')
ngiri = 3;
nxframe = 5;
nframe = ngiri*360/nxframe;
for i = 1:nframe
igrado = nxframe*i;
n = floor(igrado/360);
igrado = igrado - 360*n;
ic = find(alfa*180/pi >= igrado);
x_manov = a*cos(alfa(ic(1)));
y_manov = a*sin(alfa(ic(1)));
x_biell = c(ic(1));
y_biell = 0;
plot([0 x_manov],[0 y_manov],'r', ...
[x_manov x_biell],[y_manov y_biell],'b',...
[0 x_manov x_biell],[0 y_manov y_biell],'ok')
axis([-a-b a+b -a-b a+b])
title(sprintf('Ciclo: %d',n+1))
xlabel('[m]')
ylabel('[m]')
grid on
drawnow
end
elseif k==5
angolo=10;
ic = find(angolo/180*pi==alfa);
x_manov = a*cos(alfa(ic));
y_manov = a*sin(alfa(ic));
x_biell = c(ic(1));
y_biell = 0;
figure(10)
plot([0 x_manov],[0 y_manov],'r', ...
[x_manov x_biell],[y_manov y_biell],'b',...
[0 x_manov x_biell],[0 y_manov y_biell],'ok')
axis([-a a+b -a-b/2 a+b/2]),grid on
hold on
%plot CIR
plot([x_manov x_biell],[y_manov y_manov/x_manov*x_biell],'k:')
plot([x_biell x_biell], [0 y_manov/x_manov*x_biell],'k:')
plot(x_biell,y_manov/x_manov*x_biell,'ko')
text(x_biell*0.9,y_manov/x_manov*x_biell*1.1,'CIR')
plot([0 0],[0 y_manov/x_manov*x_biell],'k:')
plot([x_manov 0],[y_manov -tan(beta(ic))*x_biell],'k:')
plot(0,-tan(beta(ic))*x_biell,'ro')
text(-0.3,-tan(beta(ic))*x_biell*1.1,'P')
text(-0.3,0.,'O')
text(x_manov,y_manov*1.8,'A')
text(x_biell*0.95,-0.2,'B')
if abs(a+b/2)<y_manov/x_manov*x_biell
axis([-a a+b 0 y_manov/x_manov*x_biell])
end
figure(11)
subplot(2,1,1)
plot(alfa,cp), hold on
plot(alfa(ic),cp(ic),'ko'), grid
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('v_B [m/s]'))
subplot(2,1,2)
plot(alfa,-alfap*a), hold on
plot(alfa(ic),-alfap*a,'ko'), grid
xlabel('Rotazione manovella [rad]')
ylabel(sprintf('v_A [m/s]'))
ax = axis;
axis([0 2*pi ax(3) ax(4)])
grid on
% keyboard
else
stato = 0;
close all
end
end
Rispondi quotando