prandtl.m

% file PRANDTL.M : Calculation of aerodinamic characteristics of straight
% finite wings, using Prandtl - Glauert method and a Chebishev distribution
% of the points on the wing

% (C) by Luigi Tenneriello, 25/11/96
%	[e-mail];  [e-mail]
% This is a freeware program: use this program as you want,
% but don't insert it in commercial packages and don't re-sell it!

% Some important remarks:
% 1) This program was developed as an exercise in the Aerospace Engineering
%	courses of the University of Naples, Italy; I don't guarantee the
%	correctness of the results, but my teacher accepted it...; howewer,
%	it has to be (sligthly) modified to include the data of the wing you
%	want to	study in the "input section" of the program. The results
%	will be in the memory at the end of the execution.
%					   -----------------------
% 2) Note the presence of the instructions " close('all'); clear; " in the
%					   -----------------------
%	first line of the program...
% 3) Maybe the significance of the variables is obvious for an aeronautical
%	engineering student; howewer, I kept the conventions (and the
%	formula too!) from "John D. Anderson Jr.,
%	 FUNDAMENTALS OF AERODYNAMICS, Mc Graw Hill Co."
% 4) The Italian readers will be pleased to observe that:
% 	"Il programma  anche commentato in Italiano!"
% 5) Note that I wrote this program in 1 hour and an half; do you remember
%	the "Columbus egg" tale?
% 6) Thanks to the developers of Matlab!

close('all'); clear;
% INPUT SECTION
% ---------------------------------------------------------------------
n=50;			% numero di punti
			% no. of points on the wing

theta0=linspace(0,pi,n); y0=-b/2*cos(theta0);
b=5; cr=1; ct=1; S=.5*b*(cr+ct); lambda=ct/cr;
			% caratteristiche geometriche dell'ala
			% nel mio caso  un' ala trapezia
			% geometrical characteristics of the wing,
			% e.g. trapezoidal wing in my case

c=interp1([-b/2,0,b/2],[ct,cr,ct],y0)';
			% distribuzione delle corde lungo l' apertura
			% distribution of chords on the wing

AR=b^2/S,		% Allungamento Alare - Aspect Ratio

Vinf=36; rho=0.125;	% velocit / densit dell' aria
			% speed and density of the air

Clalpha=2*pi;		% in rad^-1
Clalpha=Clalpha*ones(n,1);
			% gradiente della retta di portanza di ogni profilo alare
			% lift gradient of each wing section

alpha0=0;
			% alpha_zero_lift  (in degrees)
			% di ogni profilo / of each wing section
alpha0=alpha0*ones(n,1)/180*pi;


alphar=8;		% alpha della radice (di riferimento, in gradi)
			% angle of attack of the wing root, in degrees
			% (it's the reference angle of the wing)

e=1*ones(n,1);		% svergolamento di ogni sezione, in gradi
			% twist of each wing section, in degrees


% --------------------------------------------------------------------
% END OF THE INPUT SECTION

alpha=(alphar-e)/180*pi;
			% angolo d' attacco di ogni sezione (in gradi)
			% angle of attack of each section (in degrees)

% impostazione del sistema di n-2 equazioni A*An=Anoto e risoluzione
% I build the system A*An=Anoto of n-2 equation, and my PC solve it!
Anoto=(alpha(2:n-1)-alpha0(2:n-1));
for i=2:n-1,
	i,
	for j=2:n-1,
 A(i-1,j-1)=4*b*sin((j-1)*theta0(i))/Clalpha(i)/c(i)+(j-1)*sin((j-1)*theta0(i))/sin(theta0(i));
	end;
end;
An=A\Anoto;

% calcolo di Gamma (vorticit lungo l' ala
% calculation of the vorticity Gamma on the wing
for i=2:n-1,
	Gamma(i)=2*b*Vinf*sum(An(:).*sin((1:n-2)*theta0(i))');
end;
Gamma(n)=0;Gamma(1)=0;

% calcolo altre caratteristiche aerodinamiche dell' ala
% calculation of the aerodinamical characteristics of the wing
Cl=2*Gamma./c/Vinf;
l=rho*Vinf*Gamma;	% N/m
plot(y0,c.*Cl); ylabel('c*Cl');
CL=An(1)*pi*AR,
for i=2:n-1,
	alphai(i)=sum((1:n-2)'.*An.*sin((1:n-2)'*theta0(i))./sin(theta0(i)));
end;
alphai(1)=sum((1:n-2)'.*An.*(1:n-2)');
alphai(n)=sum((1:n-2)'.*An.*(1:n-2)');
d=l.*alphai; Cdi=d./(.5*rho*Vinf^2*c);
figure; plot(y0,c.*Cdi); ylabel('c*C_d indotto');
CDi=pi*AR*sum((1:n-2)'.*(An.^2)),
L=CL*.5*rho*Vinf^2*S, Di=CDi*.5*rho*Vinf^2*S,