% návrh je provádìn pro let pøi maximální rychlosti 65 m/s s potøebným tahem 200 N

clear all
close all

R = 1.5/2;              % polomìr vrtule [m]
B = 2;                  % poèet listù vrtule [-]
V = 65;                 % rychlost [m/s]
T = 200;                % Tah [N]
%P = 15000;             % výkon
n = 2000;           	% otáèky [ot/min]
rho = 1.225;            % hustota vzduchu [kg/m^3]
cd = 0.00663;           % návrhový souèinitel odporu
cl = 0.5;               % návrhový souèinitel vztlaku
DL=cd/cl;               % pomìr souèinitele odporu a vztlaku
alfad = 0.7*pi/180;     % návrhový úhel nábìhu [°]
mi = 1.461*10^-5;       % kinematická viskozita [m^2/s]
a = 340.3;              % rychlost zvuku [m/s]

xi = 0.001:0.001:1;     % pomìrný polomìr [-]
r = xi*R;           
omega = 2*pi*n/60;      % ùhlová rychlost [rad/s]
lambda = V/(omega*R);   % inverzní rychlostní pomìr na špièce [-]
x = ((omega*xi*R)/V);   % rychlostní pomìr [-]

Tc = 2*T/((rho*V^2)*pi*R^2);     % souèinitel tahu [-]
%Pc = 2*P/(rho*V^3*pi*R^2);      % souèinitel výkonu [-]
   
f = (B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-r./R);                            % argument Prandtlovy funkce [-]
F = ((2/pi)*acos(exp(-f)));     % Prandtlova ztrátová funkce
      
G = ((F.*x.^2).*(x.^2+1).^-1);  % víøivost [-]       

% integranty i1, i2, j1, j2
i1 = ((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1-(DL)*(((omega*xi*R)/V).^-1)).*xi;

i2 = (((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1-(DL)*(((omega*xi*R)/V).^-1)).*(1*(((omega*xi*R)/V).^2+1).^-1)).*xi;

j1 = (((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1+(((omega*xi*R)/V).*DL))).*xi;

j2 = ((((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1+(((omega*xi*R)/V).*DL))).*((((omega*xi*R)/V).^2).*((((omega*xi*R)/V).^2)+1).^-1)).*xi;

% integrace pomocí Simpsonova pravidla
format long
fun1 = @(xi) ((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1-(DL)*(((omega*xi*R)/V).^-1)).*xi;
I1 = 4*quad(fun1,0,1);

fun2 = @(xi) (((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1-(DL)*(((omega*xi*R)/V).^-1)).*(1*(((omega*xi*R)/V).^2+1).^-1)).*xi;
I2 = 2*quad(fun2,0,1);

fun3 = @(xi) (((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1+(((omega*xi*R)/V).*DL))).*xi;
J1 = 4*quad(fun3,0,1);

fun4 = @(xi) ((((((2/pi)*acos(exp(-(B/2)*(sqrt(lambda^2 + 1)/lambda)*(1-xi)))).*((omega*xi*R)/V).^2).*(((omega*xi*R)/V).^2+1).^-1).*(1+(((omega*xi*R)/V).*DL))).*((((omega*xi*R)/V).^2).*((((omega*xi*R)/V).^2)+1).^-1)).*xi;
J2 = 2*quad(fun4,0,1);

zeta = (I1*(2*I2)^-1)*(1-(sqrt(1-4*Tc*I2*I1^-2)));     % rychlostní pomìr
%zeta = (J1*(2*J2)^-1)*(sqrt(1+4*Pc*J2*J1^-2)-1);      % rychlostní pomìr

Pc = J1*zeta+J2*zeta^2;                                % souèinitel výkonu
%Tc = I1*zeta-I2*zeta^2;                               % souèinite tahu

eta = Tc*Pc^-1;                                        % Úèinnost

fi = atan((lambda*xi.^-1).*(1+zeta.*0.5));      
W = sqrt(x.^2+1 - (0.5*zeta.*cos(fi)).^2)*V;                               % celková rychlost elementu listu
c = R*(4*pi*lambda*(B).^-1)*(G.*(1./((W/V)))).*(zeta*cl^-1);               % hloubka profilu elementu     

Beta = (fi + alfad)*180/pi;         % Úhel zkroucení elementu listu vrtule

Re = rho*W.*c.*(1/(mi));            % Reynoldsovo èíslo elementu listu
Ma = W.*346.3^-1;                   % Machovo èíslo elementu listu


P = 0.5*Pc*rho*V^3*pi*R^2;           % výkon
%T = Tc*rho*V^2*pi*R^2/2;            % tah

format short
display('Prùmìr vrtule '), disp(R*2)                                       
display('Poèet listù vrtule '), disp(B)                                    
display('Otáèky vrtule '), disp(n)                                         
display('Reynoldsovo èíslo v 75% délky listu '), display(Re(750))          
display('Zkroucení vrtule [°] v 75% délky listu '), display(Beta(750))     
display('Hloubka profilu [mm] v 75% délky listu '), display(1000*c(750))   
display('Machovo èíslo na špièce listu '), display(Ma(1000))               
display('Úèinnost '), disp(eta)                                            
display('Potøebný výkon [kW]'), disp(P/1000)                               

figure
plot(xi,G,'r',xi,i1,'k',xi,i2,'b',xi,j1,'g',xi,j2,'c')                     % graf závislosti G, i1, i2, j1, j2 na pomìrném polomìru xi     
grid on
xlabel('\xi [-]')
ylabel('i1, i2, j1, j2, G [-]')
axis([0,1,0,1])
legend('G', 'i1/4', 'i2/2', 'j1/4', 'j2/2')

figure                                                                     % graf závislosti úhlu zkroucení beta na pomìrném polomìru xi
plot(xi, Beta)
grid on
xlabel('\xi [-]')
ylabel('\beta [°]')
axis([0,1,0,95])

figure                                                                     % graf závislosti hloubky profilu c na pomìrném polomìru xi
plot(xi,c*1000)
grid on
xlabel('\xi [-]')
ylabel('c [mm]')
axis([0,1,0,70])

format long                       
rho2= 2770;         %Parametry materiálu Dural 2024-T3
sigmaMax =(428e6);
tmin = 0.10;        %Poèáteèní relativní tlouš�ka po celé délce listu

Ai = zeros(size(r));
dm = zeros(size(r));
dFo = zeros(size(r));
Fo = zeros(size(r));
sigma = zeros(size(r));
trel = tmin*ones(size(r));

for ii = (length(r)-1):-1:1                         %Výpoèet odstøedivé síly a kontrola dovoleného napìtí
    Ai(ii) = c(ii)^2*trel(ii);
    dm(ii) = rho2*Ai(ii)*(r(ii+1)-r(ii));
    dFo(ii) = dm(ii)*omega^2*(r(ii+1)+r(ii))/2;
    Fo(ii) = dFo(ii) + Fo(ii+1);
    sigma(ii) = 3*Fo(ii)/Ai(ii);
    if sigma(ii)>sigmaMax
        for jj = 1:10
            trel(ii) = trel(ii)*sigma(ii)/sigmaMax;
            Ai(ii) = c(ii)^2*trel(ii);
            dm(ii) = rho2*Ai(ii)*(r(ii+1)-r(ii));
            dFo(ii) = dm(ii)*omega^2*(r(ii+1)+r(ii))/2;
            Fo(ii) = dFo(ii) + Fo(ii+1);
            sigma(ii) = 3*Fo(ii)/Ai(ii);
        end
    end
end

figure                   % graf závislosti relativní tlouš�ky profilu trel na pomìrném polomìru xi
plot(xi,trel)
xlabel('\xi [-]')
ylabel('t_r_e_l [-]')
xlim([0.1 1])
grid on

figure                   % graf závislosti odstøedivé síly na pomìrném polomìru xi
plot(xi,Fo)
xlabel('\xi [-]')
ylabel('F_o [-]')
xlim([0.1 1])
grid on

figure                   % graf závislosti napìtí na pomìrném polomìru xi
plot(xi,sigma*1e-6)
xlabel('\xi [-]')
ylabel('\sigma [MPa]')
xlim([0.1 1])
grid on



