Manovellismo in Matlab

IlaCrazy
Salve a tutti,
ho un problema con Matlab: vorrei riuscire a rappresentare il meccanismo di funzionamento di un pistone sia in termini numerici che
con parte grafica. Quello che per ora sono riuscita a fare riguarda un semplice meccanismo di manovellismo ordinario centrato con sistema biella- manovella. Qualcuno riesce a darmi una mano per creare lo stesso modello per un pistone?

Grazie!!

Posto qui il manovellismo prodotto:
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
Per rispondere a questa discussione devi prima effettuare il login.