Aerodinamica

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

%% Trabajo de Aerodinámica
% Mediante el Método de los paneles se resolverá el problema de Curvatura y
% de ángulo de ataque de la linea media de un perfil. Se sacará la
% circulación de los torbellinos y con ello la velocidad de perturbación.
% Al resolver esto se podrá sacar numéricamente la sustestencación, el
% coeficiente de Momento etc...
%% Datos del Problema
fmax=0.04;
xpmax=0.25;
c=1;
Uinf=1;
alphaneg=-5*pi/180;
alpha0=0*pi/180;
alpha5=5*pi/180;
%% Calculos de los Terminos A0, A1, A2... Solo vale para este tipo de Perfil NACA
num=21;
An(1:num)=0;
barrido=10000;
k2=fmax/(1-xpmax).^2;
k1=fmax/xpmax.^2;
Int1=0;
%prompt='¿Ángulo de ataque? '
%x=double(input(prompt))
for i=(2*pi/3):1/barrido:pi
 Int1=Int1+(-0.5-cos(i))*(1/barrido);
end
Int2=0;
for i=0:1/barrido:(2*pi/3)
 Int2=Int2+(-0.5-cos(i))*(1/barrido);
end
An(1)=(1/pi)*(Int2*k2+k1*Int1);
for j=2:num;
s=0;
for i=2*pi/3:1/barrido:pi;
 s=s+(-0.5-cos(i))*cos(i*(j-1))*(1/barrido);
end
Int2=0;
for i=0:1/barrido:(2*pi)/3
 Int2=Int2+(-0.5-cos(i))*cos((j-1)*i)*(1/barrido);
end
An(j)=-(2/pi)*(s*k1+k2*Int2);
end
display(An)
%%Propiedades del Perfil
alphacl0=-An(1)-An(2)/2;
alphaid= An(1);
xca=(1/4)*(c);
alphas = [alphaid,alphaneg,alpha0,alpha5];
Cl= 2*pi*(alphas-An(1)+An(2)/2);
Cmba=-0.5*pi*(alphas-An(1)+An(2)-An(3)/2);
Xcp=0.25*c*((Cl.^(-1))*pi*(An(2)-An(3))+1);
Cmuncuarto=-0.25*pi*(An(2)+An(3));
display(Cl)
display(Cmba)
display(Xcp)
display(Cmuncuarto)
display(xca)
display(alphaid)
%% Cálculo de los paneles y sus puntos
n=100; %Ojo meter el numero de paneles y nodos que corresponda para que el punto 0.25 sea un nodo
m=n+1;
Nodos=(1:m);
Paneles=(1:n);
Xnodo=(0:1/n:c);
Znodo(1:m)=0;
for i=1:0.25*n+1
 Znodo(i)=fmax*(2*xpmax.*Xnodo(i)-Xnodo(i).^2)./(xpmax^2);
end
for i=0.25*n+2:m
 Znodo(i)=fmax*((1-2*xpmax)+2*xpmax.*Xnodo(i)-Xnodo(i).^2)/(1-xpmax)^2;
end
figure(1)
plot(Xnodo,Znodo,'b*')
hold on
plot(Xnodo,Znodo,'b')
grid
title('Linea Media de la Cuerda del Perfil')
xlabel('x/c')
ylabel('z')
%% Calculo de las Matrices del Sistema
Xci(1:n)=0;
dZ(1:n)=0;
for i=1:n
 Xci(i)= (Xnodo(i+1)+Xnodo(i)) /2;
 dZ(i)= atan((Znodo(i+1)-Znodo(i))/(Xnodo(i+1)-Xnodo(i)));
end 
B(1:n,1:n)=0;
C(1:n,1:n)=0;
for i=1:n
 for j=1:n
 rij=Xci(i)-Xnodo(j);
 rij1=Xci(i)-Xnodo(j+1);
 lj=Xnodo(j+1)-Xnodo(j);
 B(i,j)=-(1)*((lj)+(rij1)*log(abs((rij1)/(rij))))/(lj);
 if j<n
 C(i,j+1)=1+(rij*log(abs(rij1/rij))/lj);
 end
 end
end 
A(1:n,1:n)=0;
b(1:n)=0;
for i=1:n
 for j=1:n
 A(i,j)=B(i,j)+C(i,j);
 end 
end
%display(A)
%% Cálculo de los Coeficientes Aerodinámicos
%Alpha -5º
 for i=1:n
 b(i)=-2*pi*Uinf*(alphaneg-dZ(i));
 end 
 gamma=A^(-1)*b';
 gamma(m)=0;%Condicion de Kutta
 clineg(1:n)=0;
 for i=1:n
 clineg(i)=gamma(i)+gamma(i+1);
 end 
 Clneg=0;
 for i=1:n
 Clneg=(Clneg+clineg(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 cmineg(1:n)=0;
 for i=1:n
 cmineg(i)=-clineg(i)*Xci(i);
 end
 Cmneg=0;
 for i=1:n
 Cmneg=(Cmneg+cmineg(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 Xcpneg=-Cmneg/Clneg;
 
 
%Alpha 0 
 for i=1:n
 b(i)=-2*pi*Uinf*(alpha0-dZ(i));
 end 
 gamma=A^(-1)*b';
 gamma(m)=0;%Condicion de Kutta
 cli0(1:n)=0;
 for i=1:n
 cli0(i)=gamma(i)+gamma(i+1);
 end 
 Cl0=0;
 for i=1:n
 Cl0=(Cl0+cli0(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 cmi0(1:n)=0;
 for i=1:n
 cmi0(i)=-cli0(i)*Xci(i);
 end
 Cm0=0;
 for i=1:n
 Cm0=(Cm0+cmi0(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 Xcp0=-Cm0/Cl0;
%Alpha 5º
 for i=1:n
 b(i)=-2*pi*Uinf*(alpha5-dZ(i));
 end 
 gamma=A^(-1)*b';
 gamma(m)=0;%Condicion de Kutta
 cli5(1:n)=0;
 for i=1:n
 cli5(i)=gamma(i)+gamma(i+1);
 end 
 Cl5=0;
 for i=1:n
 Cl5=(Cl5+cli5(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 cmi5(1:n)=0;
 for i=1:n
 cmi5(i)=-cli5(i)*Xci(i);
 end
 Cm5=0;
 for i=1:n
 Cm5=(Cm5+cmi5(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 Xcp5=-Cm5/Cl5;
 
 
%Alpha Ideal
 for i=1:n
 b(i)=-2*pi*Uinf*(alphaid-dZ(i));
 end 
 gamma=A^(-1)*b';
 gamma(m)=0;%Condicion de Kutta
 cliid(1:n)=0;
 for i=1:n
 cliid(i)=gamma(i)+gamma(i+1);
 end 
 Clid=0;
 for i=1:n
 Clid=(Clid+cliid(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 cmiid(1:n)=0;
 for i=1:n
 cmiid(i)=-cliid(i)*Xci(i);
 end
 Cmid=0;
 for i=1:n
 Cmid=(Cmid+cmiid(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 Xcpid=-Cmid/Clid;
 
 Cl2= [Clid,Clneg, Cl0, Cl5];
 Cm2= [Cmid,Cmneg, Cm0, Cm5];
 Xcp2= [Xcpid, Xcpneg, Xcp0, Xcp5];
 display(Cl2)
 display(Cm2)
 display(Xcp2)
 alphasn = -(Cl0-alpha0*((Cl0-Clneg)/(alpha0-alphaneg)))/((Cl0-Clneg)/(alpha0-alphaneg));
 display(alphasn)
 Clalpha= (Cl0-Clneg)/(alpha0-alphaneg);
 display(Clalpha)
 
 
 
figure(2)
plot(Xci,clineg,'ro')
hold on
plot(Xci,clineg,'r')
hold on 
plot(Xci,cli0,'co')
hold on
plot(Xci,cli0,'c')
hold on 
plot(Xci,cli5,'bo')
hold on
plot(Xci,cli5,'b')
hold on
plot(Xci,cliid,'yo')
hold on
plot(Xci,cliid,'y')
grid
legend('-5º','-5º','0º','0º','5º','5º','alphaideal','alphaideal')
title('Coeficiente de Sustentación Metodo de los Paneles')
xlabel('xci/c')
ylabel('Cli') 
 
figure(3)
plot(Xci,cmineg,'r^')
hold on
plot(Xci,cmineg,'r')
hold on
plot(Xci,cmi0,'c^')
hold on
plot(Xci,cmi0,'c')
hold on
plot(Xci,cmi5,'b^')
hold on
plot(Xci,cmi5,'b')
hold on
plot(Xci,cmiid,'y^')
hold on
plot(Xci,cmiid,'y')
grid
legend('-5º','-5º','0º','0º','5º','5º','alphaideal','alphaideal')
title('Coeficiente de Momento Borde de Ataque')
xlabel('xci/c')
ylabel('Cmi')
alphas = [alphaid,alphaneg,alpha0,alpha5];
Cla = [Clid,Clneg,Cl0,Cl5];
figure(4)
plot(alphas*180/pi,Cla,'r*')
hold on
plot(alphas*180/pi,Cla,'r')
grid
title('Cl en función de alpha')
xlabel('alpha en deg')
ylabel('Cl')
figure(5)
plot(alphas*180/pi,Cm2,'r*')
hold on
plot(alphas*180/pi,Cm2,'r')
title('Cm coeficiente de momento respecto del Borde de ataque MP')
xlabel('alpha en deg')
ylabel('Cm1/4')
figure(6)
plot(alphas*180/pi,Xcp2,'r*')
title('Posición del Centro de Presiones según Alpha')
xlabel('alpha en deg')
ylabel('Xcp')
figure(7)
plot(alphas*180/pi,Cmba,'r*')
hold on
plot(alphas*180/pi,Cmba,'r')
title('Cm en el borde de ataque según TPL')
xlabel('alpha en deg')
ylabel('Cmba')
%% Cálculo de Cm1/4
 cm25id(1:n)=0;
 for i=1:n
 cmi25id(i)=-cliid(i)*(Xci(i)-0.25);
 end
 Cm25id=0;
 for i=1:n
 Cm25id=(Cm25id+cmi25id(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 cm25neg(1:n)=0;
 for i=1:n
 cmi25neg(i)=-clineg(i)*(Xci(i)-0.25);
 end
 Cm25neg=0;
 for i=1:n
 Cm25neg=(Cm25neg+cmi25neg(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 
 cm250(1:n)=0;
 for i=1:n
 cmi250(i)=-cli0(i)*(Xci(i)-0.25);
 end
 Cm250=0;
 for i=1:n
 Cm250=(Cm250+cmi250(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
 cm255(1:n)=0;
 for i=1:n
 cmi255(i)=-cli5(i)*(Xci(i)-0.25);
 end
 Cm255=0;
 for i=1:n
 Cm255=(Cm255+cmi255(i)*(Xnodo(i+1)-Xnodo(i)))/c;
 end
figure(8)
plot(Xci,cmi25neg,'r^')
hold on
plot(Xci,cmi25neg,'r')
hold on
plot(Xci,cmi250,'c^')
hold on
plot(Xci,cmi250,'c')
hold on
plot(Xci,cmi255,'b^')
hold on
plot(Xci,cmi255,'b')
hold on
plot(Xci,cmi25id,'y^')
hold on
plot(Xci,cmi25id,'y')
grid
legend('-5º','-5º','0º','0º','5º','5º','alphaideal','alphaideal')
title('Coeficiente de Momento 1/4')
xlabel('xci/c')
ylabel('Cm1/4')
CM25 = [Cm25id Cm25neg Cm250 Cm255];
display(CM25)
%% Distribución del Cl teoría potencial Linealizada 
theta(1:m+1)=0;
theta(1)=pi;
x(2:m)=Xci(1:n);
x(1)=0;
x(m+1)=1;
theta(m+1)=0;
theta(1)=pi;
for i=1:m+1
 theta(i)= acos(2*x(i)-1); 
end
cltplid(1:m+1)=0;
for i=1:m+1
 cltplid(i)=4*(alphaid-An(1))*tan(theta(i)*0.5);
 for j=2:num
 cltplid(i)= cltplid(i)+4*An(j)*sin((j-1)*theta(i)); 
 end
 
end
cltplneg(1:m+1)=0;
for i=1:m+1
 cltplneg(i)=4*(alphaneg-An(1))*tan(theta(i)*0.5);
 for j=2:num
 cltplneg(i)= cltplneg(i)+4*An(j)*sin((j-1)*theta(i)); 
 end
 
end
cltpl0(1:m+1)=0;
for i=1:m+1
 cltpl0(i)=4*(alpha0-An(1))*tan(theta(i)*0.5);
 for j=2:num
 cltpl0(i)= cltpl0(i)+4*An(j)*sin((j-1)*theta(i)); 
 end
 
end
cltpl5(1:m+1)=0;
for i=1:m+1
 cltpl5(i)=4*(alpha5-An(1))*tan(theta(i)*0.5);
 for j=2:num
 cltpl5(i)= cltpl5(i)+4*An(j)*sin((j-1)*theta(i)); 
 end
 
end
cltplid(m+1)=0;
cltplneg(m+1)=0;
cltpl0(m+1)=0;
cltpl5(m+1)=0;
figure(9)
plot(x,cltplneg,'ro') 
hold on
plot(x,cltplneg,'r') 
hold on
plot(x,cltpl0,'co') 
hold on
plot(x,cltpl0,'c') 
hold on
plot(x,cltpl5,'bo') 
hold on
plot(x,cltpl5,'b') 
hold on
plot(x,cltplid,'yo') 
hold on
plot(x,cltplid,'y') 
title('Solución por TPL')
legend('-5º','-5º','0º','0º','5º','5º','alphaideal','alphaideal')
xlabel('Xnodo')
ylabel('Cl')
axis([0 1 -3 3])
grid