clc
clear all
%% Vstupní data
D=0.105; %Vrtání [m]
r=0.060; %Poloměr kliky [m]
l=0.215; %Délka ojnice[m] 
lambda=r/l; %Ojniční poměr [-]
Sp=pi*D^2/4; %Čelní plocha pístu [m^2]
e=0.016; %Hodnota excentricity [m]
lambdae=e/l; %Excentritiký poměr [-]
ZHU=sqrt((l+r)^2-e^2); %Zdvih v HÚ [m]
ZDU=sqrt((l-r)^2-e^2); %Zdvih v DÚ [m]
Z=ZHU-ZDU; %Zdvih [m]

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]

p0=101325; %Atmosférický tlak [Pa]
n=2000; %Jmenovité otáčky motoru [min-1]
omega=(2*pi*n)/60; %Úhlová rychlost klikového hřídele[rad/s] 
alfaHU=asin(e/(l+r)); %Poloha kliky při HÚ [rad]
angleHU=asin(e/(l+r))*(180/pi); %Poloha kliky při HÚ [°]
alfaDU=asin(e/(l-r))+pi; %Poloha kliky při DÚ [rad]
angleDU=(asin(e/(l-r))+pi)*(180/pi); %Poloha kliky při DÚ [°]
%% Načtení dat 
%Tlaky
% p_input=xlsread('\data\tlaky.xlsx','PosunuteTlaky','R2:R721');
% p = interp1(1:720,p_input,1:0.01:720);
 p=xlsread('\data\tlaky.xlsx','PosunuteTlaky','R2:R721');
%Úhel natočení kliky
% angle=1:0.01:720; %[°]
angle=1:1:720; %[°]
alfa=(angle*pi)/180; %[rad]
%% Výpočty ztrátového výkonu
for i=1:length(alfa)
        s(i)=sqrt((l+r)^2-e^2)-r*(1/lambda+cos(alfa(i))+lambdae*sin(alfa(i))-lambda/4*(1-cos(2*alfa(i)))); %Dráha pístu [m]    
        s1(i)=-r*cos(alfa(i))-r*lambdae*sin(alfa(i)); %1.harmonická složka [m]
        s2(i)=-r*lambda/4*cos(2*alfa(i)); %2.harmonická složka [m]
        
        v(i)=r*omega*(sin(alfa(i))-lambdae*cos(alfa(i))+lambda/2*sin(2*alfa(i))); %Rychlost pístu [m/s]  
        v1(i)=r*omega*(sin(alfa(i))-lambdae*cos(alfa(i))); %1.harmonická složka [m/s]
        v2(i)=r*omega*lambda/2*sin(2*alfa(i)); %2.harmonická složka [m/s]
                              
        a(i)=r*omega^2*(cos(alfa(i))+lambdae*sin(alfa(i))+lambda*cos(2*alfa(i))); %Zrychlení pístu [m/s2] 
        a1(i)=r*omega^2*(cos(alfa(i))+lambdae*sin(alfa(i))); %1.harmonická složka [m/s2]
        a2(i)=r*omega^2*lambda*cos(2*alfa(i)); %2.harmonická složka [m/s2]
           
        Fp=Sp*(p*100000-p0)'; %Síly od tlaku plynů [N]
        Fs(i)=-mp*a(i); %Setrvačné síly pohybujících se částí [N]
        Fc(i)=Fp(i)+Fs(i); %Celková síla [N]
end
for i=1:720
    cs(i)=2*Z*n/60; %Střední pístová rychlost [m/s2]
end

cs_plus=mean(v(3:185)); %Střední pístová rychlost v horním zdvihu [m/s2]
cs_minus=mean(v(185:363)); %Střední pístová rychlost v dolním zdvihu [m/s2]

Fnc=0;
for i=1:length(alfa)
    beta(i)=asin(lambda*sin(alfa(i))-lambdae); %Úhel odklonu ojnice [rad]
    Fn(i)=Fc(i)*tan(beta(i)); %Normálová síla [N]
end
Fnc=mean(Fn); %Střední hodnota normálové síly [N]
disp(['Střední hodnota síly mezi pístem a vložkou válce e = 16mm (2000 min-1) = ',num2str(Fnc), ' N'])
Pzc=0;
for i=1:length(alfa)
    Pz(i)=abs(Fn(i))*abs(v(i))*0.05; %Ztrátový výkon [W]
end
Pzc=mean(Pz)*2; %Celková hodnota ztrátového výkonu(* počet válců)[W]
disp(['Ztrátový výkon dvouválce e = 16mm (2000 min-1) = ',num2str(Pzc), ' W'])
%% Výpočet točivého a klopného momentu
Fsr=-mro*r*omega^2; %Setrvačná síla rotujících částí [N]
Mtc=0;
for i=1:length(alfa)
    Fo(i)=Fc(i)/cos(beta(i)); %Ojniční síla [N]
    Ftan(i)=(Fo(i))*sin(alfa(i)+beta(i)); %Tangenciální síla [N] 
    Fr(i)=(Fo(i))*cos(alfa(i)+beta(i)); %Radiální síla [N]
    Fcr(i)=Fr(i)+Fsr; %Celková radiální síla [N]
    Mt(i)=Ftan(i)*r; %Točivý moment [Nm]
end
Mtc=mean(Mt)*2; %Celková hodnota točivého momentu dvouválce [Nm]
disp(['Točivý moment dvouválce e = 16mm (2000 min-1) = ',num2str(Mtc), ' Nm'])

Mklc=0;   
for i=1:length(alfa)
    b(i)=r*(sin(alfa(i)+beta(i))/sin(beta(i))); %Vzdálenost od osy pístního čepu k ose kliky [m]
    Mkl(i)=Fn(i)*b(i); %Klopný moment [Nm]
end
Mklc=mean(Mkl)*2; %Celkový klopný moment [Nm]

