Importing different plots onto one plot
4 views (last 30 days)
Show older comments
I have a matlab code that prompts the user to enter specific details about an airfoil. then the user is able to plot different aspects of the airfoil. what i need to do is combine the plots from these two separate runs of the program. i have up loaded the entire code. I'm trying to explain it the best I can. basically, when you run the code for one type of airfoil, you get a certain result. when u run the code again for a different type of airfoil, new results are found, but the results from the initial airfoil are then gone. Is there a way to take the results from separate runs of the program and combine them into one plot? Thanks
clear all
clc
disp('Running MAE 424 vortex panel method code.');
disp(' ');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part I: user input (airfoil geometry, # panels, min/max angles of attack)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pmxx = input('Enter desired NACA 4-digit designation, e.g. 0012: ');
disp(' ');
Np = input('Enter number of vortex panels (even integer, ideally 40-100): ');
disp(' ');
alpha_min = input('Enter the minimum angle of attack to test (in degrees): ');
disp(' ');
alpha_max = input('Enter the maximum angle of attack to test (in degrees): ');
disp(' ');
kmax = input('Enter the max number of angles to test over this range (integer): ');
disp(' ');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part II: creation of NACA 4-digit airfoil panely geometry
%
% Calculates the coordinates of the airfoil panel boundary points
% (xfoil,yfoil), and coordinates of the mean camber line points (xc,yc)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
disp('Creating (x,y) coordinates of airfoil boundary points, camber line...');
disp(' ');
pm = fix(pmxx/100);
xx = pmxx-pm*100;
p = fix(pm/10);
m = pm-p*10;
p = p/100;
m = m/10;
if m~=0
fact1 = p/m^2; fact2 = 2*m; fact3 = p/(1-m)^2;
else
fact1 = 0; fact2 = 19; fact3 = 0;
end
acoef = xx/20*[0.2969 -0.1260 -0.3516 0.2843 -0.1015];
%
% Compute thickness distribution
% vertices are clustered near the leading (x/c=0) and
% trailing (x/c=1) edges
%
beta = 0:2*pi/Np:2*pi; xc = 0.5*(1+cos(beta));
yt = acoef(1)*sqrt(xc)+xc.*(acoef(2)+xc.*(acoef(3)+...
xc.*(acoef(4)+acoef(5)*xc)));
yt(Np/2+2:Np+1) = -yt(Np/2:-1:1);
%
% Add the mean camber line
%
for i=1:Np+1
if xc(i)<m
yc(i) = fact1*xc(i)*(fact2-xc(i));
tand = fact1*(fact2-2*xc(i));
else
yc(i) = fact3*(1-fact2+xc(i)*(fact2-xc(i)));
tand = fact3*(fact2-2*xc(i));
end
cosd = 1/sqrt(1+tand*tand); sind = tand*cosd;
xfoil(i) = xc(i)+yt(i)*sind; yfoil(i) = yc(i)-yt(i)*cosd;
end
% Plot airfoil panel geometry to check
figure(1)
plot(xfoil,yfoil,'r');
daspect([1 1 1]);
hold on
plot(xc,yc,'b');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part III: calculations using vortex panel method
%
% Calculates: gamma - array of circulation densities;
% xcen - ordinates of panel centers;
% cp - pressure coefficient;
% clift - lift coefficient;
% cdrag - drag coefficient (aerodynamic);
% ista - starting panel i index;
% iend - ending panel i index;
% xcp - center of pressure;
% cmle - moment coeff. about leading edge;
% cmqc - moment coeff. about 1/4 chord.
%
% All notations follow Kuethe & Chow (1998)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
disp('Running vortex panel method...');
disp(' ');
Np1 = Np+1;
alpha_min = alpha_min*pi/180;
alpha_max = alpha_max*pi/180;
alpha = linspace(alpha_min,alpha_max,kmax)'; % Angle of attack vector
alpha_deg = alpha*180/pi;
% Initialize column vectors of calculated quantities, fill with 0's
clift = zeros(kmax,1);
cdrag = zeros(kmax,1);
xcp = zeros(kmax,1);
cmle = zeros(kmax,1);
cmqc = zeros(kmax,1);
cl_exp = [.00395, .204, .308, .613, .917, 1.2065,1.4591];
a_exp = [.138,2.0496,2.9303,6.0206,8.8127,12.1,15];
for k = 1:kmax % Overall 'for' loop for varying angle of attack
%
% Compute geometrical characteristics
%
for i=1:Np,
ip1 = i+1;
x(i) = 0.5*(xfoil(i)+xfoil(ip1));
y(i) = 0.5*(yfoil(i)+yfoil(ip1));
s(i) = sqrt((xfoil(i)-xfoil(ip1))^2+...
(yfoil(i)-yfoil(ip1))^2);
theta(i) = atan2((yfoil(ip1)-yfoil(i)),...
(xfoil(ip1)-xfoil(i)));
sine(i) = sin(theta(i)); cosine(i) = cos(theta(i));
rhs(i) = sin(theta(i)-alpha(k));
end
%
% Assemble the system matrix
%
for i=1:Np, for j=1:Np
if (i==j)
cn1(i,j)=-1;cn2(i,j)=1;ct1(i,j)=0.5*pi;ct2(i,j)=0.5*pi;
else
A =-(x(i)-xfoil(j))*cosine(j)-(y(i)-yfoil(j))*sine(j);
B = (x(i)-xfoil(j))^2+(y(i)-yfoil(j))^2;
C = sin(theta(i)-theta(j));
D = cos(theta(i)-theta(j));
E = (x(i)-xfoil(j))*sine(j)-(y(i)-yfoil(j))*cosine(j);
F = log(1+s(j)*(s(j)+2*A)/B);
G = atan2(E*s(j),B+A*s(j));
P = (x(i)-xfoil(j))*sin(theta(i)-2*theta(j))+...
(y(i)-yfoil(j))*cos(theta(i)-2*theta(j));
Q = (x(i)-xfoil(j))*cos(theta(i)-2*theta(j))-...
(y(i)-yfoil(j))*sin(theta(i)-2*theta(j));
cn2(i,j) = D+0.5*Q*F/s(j)-(A*C+D*E)*G/s(j);
cn1(i,j) = 0.5*D*F+C*G-cn2(i,j);
ct2(i,j) = C+0.5*P*F/s(j)+(A*D-C*E)*G/s(j);
ct1(i,j) = 0.5*C*F-D*G-ct2(i,j);
end
end, end
%
%
%
for i=1:Np
An(i,1) = cn1(i,1); An(i,Np1) = cn2(i,Np);
At(i,1) = ct1(i,1); At(i,Np1) = ct2(i,Np);
for j=2:Np
An(i,j)=cn1(i,j)+cn2(i,j-1); At(i,j)=ct1(i,j)+ct2(i,j-1);
end
end
An(Np1,1) = 1; An(Np1,Np1) = 1;
for j=2:Np, An(Np1,j) = 0; end
rhs(Np1) = 0;
%
% Solve for vortex panel strengths
%
gamma = An\rhs';
%
xcen = x;
i=1; while xcen(i+1)>0.99, i=i+1; end
ista = i;
i=Np; while xcen(i-1)>0.99, i=i-1; end
iend = i;
%
% Compute Cp, Clift, and Cdrag
%
cx = 0; cy = 0; cmle(k) = 0; cmqc(k) = 0;
for i=ista:iend
v(i) = cos(theta(i)-alpha(k));
for j=1:Np1
v(i) = v(i)+At(i,j)*gamma(j);
end
cp(i) = 1-v(i)^2;
cx = cx+cp(i)*s(i)*sine(i);
cy = cy-cp(i)*s(i)*cosine(i);
cmle(k) = cmle(k)+cp(i)*s(i)*cosine(i)*xcen(i);
cmqc(k) = cmqc(k)+cp(i)*s(i)*cosine(i)*(xcen(i)-0.25);
end
clift(k) = cy*cos(alpha(k))-cx*sin(alpha(k));
cdrag(k) = cy*sin(alpha(k))+cx*cos(alpha(k));
xcp(k) = -cmle(k)/clift(k);
end % for k = 1:
% Basic plot of lift coefficient vs. angle of attack
figure(2)
hold on
plot(alpha_deg,clift,'r')
xlabel('Alpha (deg)')
figure(2)
hold on
plot(alpha_deg,cdrag,'g')
figure(2)
hold on
plot(alpha_deg,cmqc,'b')
figure(2)
hold on
plot(alpha_deg,xcp,'m')
figure(2)
hold on
plot(a_exp,cl_exp,'k')
title('NACA 0012 - C_{drag},C_{lift},C_{m,1/4},X_{CP},C_{lift,exp} vs Alpha')
ylabel('C_{lift},C_{drag},C_{m,1/4},X_{cp},C_{lift,exp}')
legend('C_{lift}','C_{drag}','C_{m/1/4}','X_{cp}','C_{lift,exp}')
1 Comment
per isakson
on 7 May 2016
There are many contributions in the File exchange, which address this problem, e.g. see Panel by Ben Mitch
Answers (1)
See Also
Categories
Find more on Coordinate Systems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!