clc
clear all
%% Vstupní data
r=0.060; %Poloměr kliky [m]
l=0.215; %Délka ojnice[m] 
lambda=r/l; %Ojniční poměr [-]

Dr=0.08; %Průměr redukovaný [m]
Dh=0.08; %Vnější průměr hlavního čepu [m]
Do=0.066; %Vnější průměr ojničního čepu [m]
lh=0.044; %Délka hlavního čepu [m]
lo=0.04; %Délka ojničního čepu [m]
lram=0.026; %Délka jednoho ramena zalomení [m]
h=0.106; %Šířka zalomení kliky [m] 
l_zalomeni=0.136; %Délka zalomení [m]

G=81000*1000000; %Modul pružnosti ve smyku [MPa;N*m^-2]

mpo=0.899; %Hmotnost posuvné části ojnice [kg]
m1=1.364; %Hmotnost pístu [kg]
m2=0.6138; %Hmotnost pístního čepu [kg]
m3=2*0.00483; %Hmotnost pojistných kroužků [kg]
m4=0.0233; %Hmotnost 1 těsnícího kroužku [kg]
m5=0.0320; %Hmotnost 2 těsnícího kroužku [kg]
m6=0.03; %Hmotnost stíracího kroužku [kg]
mp=mpo+m1+m2+m3+m4+m5+m6; %Hmotnost posuvných částí [kg] 
mro=1.662; %Hmotnost rotující části ojnice [kg]

Jzal=0.0453; %Moment setrvačnosti jednoho zalomení [kg*m^2]
Jr=0.02; %Moment setrvačnosti řemenice [kg*m^2]
Jhr=0.0027; %Moment setrvačnosti hřídele od hlavního čepu po řemenici [kg*m^2]
Js=1; %Moment setrvačnosti setrvačníku [kg*m^2]
Jhs=0.0052; %Moment setrvačnosti hřídele od hlavního čepu po setrvačník [kg*m^2]

%% Redukce hmotností
Jro=mro*r^2; %Moment setrvačnosti rotující části ojnice [kg*m^2]
Jp=mp*(1/2+lambda^2/8)*r^2; %Moment setrvačnosti posuvných hmot [kg*m^2]

Ji=Jzal+Jro+Jp; %Celkový moment setrvačnosti jednoho zalomení [kg*m^2]
disp(['Celkový moment setrvačnosti zalomení = ',num2str(Ji), ' kg*m^2'])

J1=Jr+Jhr; %Moment setrvačnosti na straně řemenice [kg*m^2]
disp(['Moment setrvačnosti u řemenice = ',num2str(J1), ' kg*m^2'])

J4=Js+Jhs; %Moment setrvačnosti na straně setrvačníku [kg*m^2]
disp(['Moment setrvačnosti u setrvačníku = ',num2str(J4), ' kg*m^2'])

%% Redukce délek
% Redukce délky jednoho zalomení [m]
lred_z=Dr^4*(((lh+(0.4*Dh))/Dh^4)+((lo+(0.4*Do))/Do^4)+((r-(0.2*(Dh+Do)))/lram*h^3));
disp(['Redukovaná délka zalomení = ',num2str(lred_z), ' m'])

% Redukce délky u řemenice [m]
lred_r_1=0.05*(Dr^4/(0.045^4*(1-((0.022/0.045)^4))));
lred_r_2=0.046*(Dr^4/0.045^4);
lred_r_3=0.025*(Dr/0.0421)^4;
lred_r_4=((0.008+(0.093*0.045))*(Dr^4/0.045^4))+(0.022-(0.093*0.045)*(Dr^4/0.06^4));
lred_r=lred_r_1+lred_r_2+lred_r_3+lred_r_4+lred_z/2;
disp(['Redukovaná délka u řemenice = ',num2str(lred_r), ' m'])

% Redukce délky u setrvačníku [m]
lred_s_1=(0.022+(0.07*0.08)*(Dr^4/0.08^4))+(0.020-(0.07*0.110)*(Dr^4/0.110^4));
lred_s_2=0.015*(Dr^4/(0.110^4*(1-(0.052/0.110)^4)));
lred_s=lred_s_1+lred_s_2+lred_z/2;
disp(['Redukovaná délka u setrvačníku = ',num2str(lred_s), ' m'])

%% Výpočet tuhostí
Ip=pi*Dr^4/32; %Polární kvadratický moment [m^4]
k1=G*Ip/lred_r; %Torzní tuhost u řemenice [Nmrad^-1]
k2=G*Ip/lred_z; %Torzní tuhost zalomení [Nmrad^-1]
k3=G*Ip/lred_s; %Torzní tuhost u setrvačníku [Nmrad^-1]
disp(['Torzní tuhost u řemenice = ',num2str(k1), ' Nmrad^-1'])
disp(['Torzní tuhost mezi zalomeními = ',num2str(k2), ' Nmrad^-1'])
disp(['Torzní tuhost u setrvačníku = ',num2str(k3), ' Nmrad^-1'])

%% Výpočet vlastních frekvencí
M=[J1 0 0 0; 0 Ji 0 0; 0 0 Ji 0; 0 0 0 J4]; %Matice hmotností
K=[k1 -k1 0 0; -k1 k1+k2 -k2 0; 0 -k2 k2+k3 -k3; 0 0 -k3 k3]; %Matice tuhostí
As=eig(inv(M)*K);
for i=1:3
OMEGA(i)=sqrt(abs(As(i)))/(2*pi); %/2*pi je prevod rad na Hz
end
disp(['První vlastní frekvence = ',num2str(OMEGA(3)), ' Hz'])
disp(['Druhá vlastní frekvence = ',num2str(OMEGA(2)), ' Hz'])

%% Výpočet vynuceného kmitání
kappa=0.5:0.5:12;%Řád harmonické složky [-]
N(1)=OMEGA(3);
N(2)=OMEGA(2);
for i=1:2
    nr(:,i)=60*N(i)./kappa; %Rezonanční otáčky jednouzlového a dvojuzlového kmitání %Řád harmonické složky [min-1]
end


