lifting_line

%%%%%%%%%%%%% TAREFA 2.2 - LIFTING LINE %%%%%%%%%%%%%

%           Augusto Morcelli e Jonas Muller.


Npanels = 90; %número de paineis
ct = 100; %corda na ponta (tip chord)[mm]
cr = 150; %corda do meio[mm]
b = 1000; %envergadura[mm]
S = ((cr+ct)*b/2)*10^(-6); %área alar [m^2]
AR = b^2/S; %Aspect Ratio (alongamento)
y = -b/2:b/Npanels:b/2;
rho_o = 1.1839; %densidade do escoamento a 25ºC[Kg/m^3]
V_o = 3.5; %Vel. inicial do escoamento [m/s]
q_o = 0.5*rho_o*(V_o^2); %pressão dinâmica
n = 1; %fator de carga
m = 1000; %massa [kg]
g = 9.81; %aceleração da gravidade [m/s^2]
W = m*g; %Peso [N]
L_o = n*W; %sustentação
circ_o = 4*L_o/(rho_o*V_o*b*pi); %chute inicial
enflex = atan2((cr/2-ct/2),(y(Npanels/2+1)-y(1))); %angulo de enflexamento
enflex_dg = enflex*180/pi; %angulo de enflexamento em graus


% enflex_dg = 180; %angulo de enflexamento em graus
% enflex = enflex_dg*pi/180; %angulo de enflexamento



circ = circ_o*sqrt((1-(2.*y./b).^2));

plot(y,circ);

% x = 2*y/b*ct/2+cr/2


for i = 1:Npanels+1

    if i == Npanels+1
        break
    end

    dcirc_dy(i) = (circ(i+1)-circ(i))./(y(i+1)-y(i));
    dy(i) = y(i+1)-y(i);

end


for i = 1:Npanels

    for j = 1:Npanels

        alpha_indaux(i,j) =(1/(4*pi*V_o)).*dcirc_dy(j).*dy(i)./(y(i) - y(j));

        if j == i

                if j == 1
                    alpha_indaux(i,j) = (1/(4*pi*V_o)).*dcirc_dy(j+1).*dy(i)./(y(i) - y(j+1));
                end

                if j == Npanels
                    alpha_indaux(i,j) = (1/(4*pi*V_o)).*dcirc_dy(j-1).*dy(i)./(y(i) - y(j-1));
                end

                if (j > 1) && (j < Npanels)

                    alpha_indaux(i,j) = (alpha_indaux(i,j-1)+alpha_indaux(i,j+1))./2;

                end

        end
    end
end

alpha_indaux_dg = alpha_indaux*180/pi;

alpha_ind = sum(alpha_indaux,2);
alpha_ind_dg = alpha_ind*180/pi;


for i = 1:Npanels
    CL_aux(i) = (2/(V_o*S))*circ(i)*dy(i);
    CDi_aux(i) = (2/(V_o*S))*circ(i)*alpha_ind(i).'*dy(i); %.' faz a transposta do vetor para transformar de vetor coluna para vetor linha.
    CDi2_aux(i) = (circ_o*b*pi/(2*V_o*S))*alpha_ind(i).';
    alpha_ind2(i) = CL_aux(i)/(pi*AR);
    L_aux(i) = rho_o*V_o*circ(i)*dy(i); %<<< distribuição de sustentação
    Di_aux(i) = rho_o*V_o*circ(i)*alpha_ind(i).'*dy(i);
    Di2_aux(i) = L_aux(i).'*alpha_ind(i);
end

CL = sum(CL_aux); %<<<<< estimativa de CL
CDi = sum(CDi_aux); %<<<<< estimativa de CDi

CDi2 = sum(CDi2_aux); %para conferir [XXX]
CDi3 = (CL^2)/(pi*AR); %para conferir
circ_o2 = 2*CL*V_o*S/(b*pi); %para conferir
L = sum(L_aux);%para conferir
Di = sum(Di_aux);
Di2 = sum(Di2_aux);%para conferir
CL2 = L/(q_o*S); %para conferir


%%%% wingplot %%%%

coordinates_up(:,1) = linspace(-b/2, b/2,Npanels+1);
coordinates_up(1:(Npanels/2 +1),2) = linspace(ct/2, cr/2,(Npanels/2)+1);
coordinates_up((Npanels/2 +1):(Npanels+1),2) = linspace(cr/2,ct/2,(Npanels/2)+1);


coordinates_down(:,1) = coordinates_up(:,1);
coordinates_down(:,2) = -coordinates_up(:,2);

for i = 1:Npanels+1

    if i == Npanels+1
        break
    end

    faixa(i,1) = (y(i+1)+y(i))/2;
    faixa(i,2) = ((coordinates_up(i,2)+coordinates_down(i,2)));
   % dy(i) = y(i+1)-y(i);

end

zaxis = zeros(length(faixa),1);

figure(2)
plot(coordinates_up(:,1),coordinates_up(:,2),'bo')
hold on;
plot(coordinates_down(:,1),coordinates_down(:,2),'ro')
grid on; grid minor;

quiver3(faixa(:,1),zaxis, faixa(:,2),zaxis,zaxis, L_aux.',0.1); %plot da distribuição de sustentação