## how do i rotate an airfoil points? and plot.

### Jorge Cabrera (view profile)

on 18 Apr 2019
Latest activity Edited by David Wilson

on 19 Apr 2019

### David Wilson (view profile)

How do I rotate an airfoil points? I have all the new chord lengths and the angles in which are needed to rotate. I need to plot all these airfoils 16 times together in respect to their new angle and chord length. I know I need a sort of for loop but evetrytime I plot, the figure is funky looking.
NEW Chord lengths and twist angle, R vaules
Chord length= [ 0 0.0082328 0.0164657 0.0246985 0.0329313 0.0411642 0.049397 0.0576298 0.0658626 0.07409548 0.08232831 0.09056114 0.09879397 0.1070268 0.11525963 0.12349246 ]
r =[ 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 0.33 0.36 0.39 0.42 0.45]
twist = [1.0471976 0.7561128 0.5465782 0.4134997 0.3278728 0.2699279 0.2286826 0.1980428 0.1744757 0.15582879 0.14072889 0.12826373 0.1178059 0.10891077 0.10125497 0.09459804]
BELOW
THIS IS THE CODE I STATRED WITH FOR NACA
THE PLOT SHOULD LOOK SOMETHING LIKE THIS BELOW
c1= 1;
%c1= [0,0.008232831,0.016465661,0.024698492,0.032931323,0.041164153,0.049396984,0.057629815,0.065862645,0.074095476,0.082328307,0.090561137,0.098793968,0.107026798,0.115259629,0.12349246];
angle =[1.047197551 0.756112778 0.546578176 0.413499657 0.327872784 0.269927858 0.228682627 0.198042808 0.174475668 0.155828787 0.140728889 0.128263734 0.117805904 0.108910765 0.101254972 0.094598036];
s=num2str(5410);
%4 digits
NACA=s;
% m represents first digit out of the scalar
m =str2double(s(1))/100;
%p pulls the second digit out of the scalar
p=str2double(s(2))/10;
%t pulls and represents the third and fourth digit out of the scalar
t=str2double(s(3:4))/100;
% define x
r=linspace(0, c1, 250);
%yt fuction including the coeffiencients
yt =5*t*c1*(.2969*(sqrt(r/c1))+-.1260*(r/c1)+-.3516*(r/c1).^2+.2843*(r/c1).^3+-.1015*(r/c1).^4);
%yt2 =5*t*c2*(.2969*(sqrt(x/c2))+-.1260*(x/c2)+-.3516*(x/c2).^2+.2843*(x/c2).^3+-.1015*(x/c2).^4);
for k = 1:length(r)
if r(k) <= p*c1
yc(k)=m*(r(k)/p^2)*(2*p-(r(k)/c1));
dx(k)=(2*m)/p^2*(p-(r(k)/c1));
elseif r(k) > p*c1
yc(k)=m*((c1-r(k))/(1-p)^2)*(1+(r(k)/c1)-(2*p));
dx(k)=((2*m)/(1-p)^2)*(p-(r(k)/c1));
end
%upper and lower limits of the airfoil (xu,yu) ; (xl,yl)
angle=atan(dx(k));
xu(k)=r(k)-yt(k)*-sin(angle);
yu(k)=yc(k)+yt(k)*cos(angle);
xl(k)=r(k)+yt(k)*sin(angle);
yl(k)=yc(k)-yt(k)*cos(angle);
end
%plot of airfoil
plot(xu,yu,xlabel('x'),ylabel('y'));
title( ['NACA Airfoil ', NACA]);
hold on
plot(xl,yl,'r')
plot(r,yc,'g')
grid on
axis equal
%--------------------------------------------------------------------------------------%
% Given data vectors X and Y.
% Want to rotate the data through angle "ang" about rotation center Xc, Yc
X = [xl,xu];
Y = [yl,yu];
ang = [254]; % 0.756112778 0.546578176 0.413499657 0.327872784 0.269927858 0.228682627 0.198042808 0.174475668 0.155828787 0.140728889 0.128263734 0.117805904 0.108910765 0.101254972 0.094598036];
% Specify the coordinates of the center of rotation
Xc = 0.25 ; % Rotate about the 1/4 chord point
Yc = 0 ;
% The data is roated in a three-step process
% Step 1) Shift the data to the rotation center
Xs = X - Xc; % shifted data
Ys = Y - Yc;
% Step 2) Rotate the data
Xsr = Xs.*cos(ang) + Ys.*sin(ang); % shifted and rotated data
Ysr = -Xs*sin(ang) + Ys*cos(ang); %
% Step 3) Un-shift the data (back to the original coordinate system)
Xr = Xsr + Xc; % Rotated data
Yr = Ysr + Yc;
hold on

### David Wilson (view profile)

on 18 Apr 2019
Edited by David Wilson

### David Wilson (view profile)

on 18 Apr 2019

Does this look better?
Some (small) changes to your code. I swapped your rotation direction to better follow your example plot. Probably not a good idea to call a variable angle.
c1= 1;
%c1= [0,0.008232831,0.016465661,0.024698492,0.032931323,0.041164153,0.049396984,0.057629815,0.065862645,0.074095476,0.082328307,0.090561137,0.098793968,0.107026798,0.115259629,0.12349246];
s=num2str(5410);
%4 digits
NACA=s;
% m represents first digit out of the scalar
m =str2double(s(1))/100;
%p pulls the second digit out of the scalar
p=str2double(s(2))/10;
%t pulls and represents the third and fourth digit out of the scalar
t=str2double(s(3:4))/100;
% define x
r=linspace(0, c1, 250);
%yt fuction including the coeffiencients
yt =5*t*c1*(.2969*(sqrt(r/c1))+-.1260*(r/c1)+-.3516*(r/c1).^2+.2843*(r/c1).^3+-.1015*(r/c1).^4);
%yt2 =5*t*c2*(.2969*(sqrt(x/c2))+-.1260*(x/c2)+-.3516*(x/c2).^2+.2843*(x/c2).^3+-.1015*(x/c2).^4);
for k = 1:length(r)
if r(k) <= p*c1
yc(k)=m*(r(k)/p^2)*(2*p-(r(k)/c1));
dx(k)=(2*m)/p^2*(p-(r(k)/c1));
elseif r(k) > p*c1
yc(k)=m*((c1-r(k))/(1-p)^2)*(1+(r(k)/c1)-(2*p));
dx(k)=((2*m)/(1-p)^2)*(p-(r(k)/c1));
end
%upper and lower limits of the airfoil (xu,yu) ; (xl,yl)
angle=atan(dx(k));
xu(k)=r(k)-yt(k)*-sin(angle);
yu(k)=yc(k)+yt(k)*cos(angle);
xl(k)=r(k)+yt(k)*sin(angle);
yl(k)=yc(k)-yt(k)*cos(angle);
end
%% plot of airfoil
clf;
plot(xu,yu,xl,yl,r,yc);
xlabel('x'),ylabel('y')
title( ['NACA Airfoil ', NACA]);
grid on; axis equal
%% --------------------------------------------------------------------------------------%
% Given data vectors X and Y.
% Want to rotate the data through angle "ang" about rotation center Xc, Yc
X = [xl,fliplr(xu)]; % watch this
Y = [yl,fliplr(yu)];
% Specify the coordinates of the center of rotation
Xc = 0.25 ; % Rotate about the 1/4 chord point
Yc = 0 ;
% The data is roated in a three-step process
% Step 1) Shift the data to the rotation center
Xs = X - Xc; % shifted data
Ys = Y - Yc;
% Step 2) Rotate the data
angles =[1.047197551 0.756112778 0.546578176 0.413499657 0.327872784 ...
0.269927858 0.228682627 0.198042808 0.174475668 0.155828787 0.140728889 ...
0.128263734 0.117805904 0.108910765 0.101254972 0.094598036];
hold on
for i=1:length(angles)
ang = -angles(i); % in radians ??
Xsr = Xs.*cos(ang) + Ys.*sin(ang); % shifted and rotated data
Ysr = -Xs*sin(ang) + Ys*cos(ang); %
% Step 3) Un-shift the data (back to the original coordinate system)
Xr = Xsr + Xc; % Rotated data
Yr = Ysr + Yc;
plot(Xr, Yr)
end % for
hold off

Jorge Cabrera

### Jorge Cabrera (view profile)

on 18 Apr 2019
Thank you David! plot works fine. But how would i get it to plot with new chord lenght values?
c2= [0,0.008232831,0.016465661,0.024698492,0.032931323...
0.041164153,0.049396984,0.057629815,0.065862645...
0.074095476,0.082328307,0.090561137,0.098793968,0.107026798...
0.115259629,0.12349246];
to plot with the same rotated angles
David Wilson

### David Wilson (view profile)

on 19 Apr 2019
Like this?
I've rotated each of your 16 profiles by the corresponding 16 angles. I assume this is what you want. By the way, you do have a chord length of zero! Seems a bit strange.
I've rotated the aerofoils about 1/4 chord length of the actual chord. Looking at your plot above, I see that you have probably rotated about the point 1/4 max chord length. It's an easy change to make if that's what you want.
I've also lightly coloured the foils to make them look nicer.
Code is now as below. I've added a structure to contain each of the profiles. This could be in a rectangular matrix, but you might in the future want to add other fields. (BTW, I typically use English spelling, simply to piss some people off.)
%% Rotate NACA aerofoils with adjusting chord lengths
s=num2str(5410); NACA=s; %4 digits
m =str2double(s(1))/100; % m represents first digit out of the scalar
p=str2double(s(2))/10; %p pulls the second digit out of the scalar
%t pulls and represents the third and fourth digit out of the scalar
t=str2double(s(3:4))/100;
% Chord length, but dropping first (zero) entry
c2= [0,0.008232831,0.016465661,0.024698492,0.032931323...
0.041164153,0.049396984,0.057629815,0.065862645...
0.074095476,0.082328307,0.090561137,0.098793968,0.107026798...
0.115259629,0.12349246];
%c2(1) = []; % drop the zero chord length !!
aerofoil = struct('chord',[],'x',[],'y',[],'r',[],'yc',[]); % place holder
for i=1:length(c2)
c1 = c2(i);
r=linspace(0, c1, 250)'; % define x
%yt function including the coefficients
yt =5*t*c1*(.2969*(sqrt(r/c1))+-.1260*(r/c1)+-.3516*(r/c1).^2+.2843*(r/c1).^3+-.1015*(r/c1).^4);
%yt2 =5*t*c2*(.2969*(sqrt(x/c2))+-.1260*(x/c2)+-.3516*(x/c2).^2+.2843*(x/c2).^3+-.1015*(x/c2).^4);
for k = 1:length(r)
if r(k) <= p*c1
yc(k)=m*(r(k)/p^2)*(2*p-(r(k)/c1));
dx(k)=(2*m)/p^2*(p-(r(k)/c1));
elseif r(k) > p*c1
yc(k)=m*((c1-r(k))/(1-p)^2)*(1+(r(k)/c1)-(2*p));
dx(k)=((2*m)/(1-p)^2)*(p-(r(k)/c1));
end
%upper and lower limits of the airfoil (xu,yu) ; (xl,yl)
angle=atan(dx(k));
xu(k)=r(k)-yt(k)*-sin(angle);
yu(k)=yc(k)+yt(k)*cos(angle);
xl(k)=r(k)+yt(k)*sin(angle);
yl(k)=yc(k)-yt(k)*cos(angle);
end % for k
aerofoil(i).chord = c1;
aerofoil(i).x = [xl,fliplr(xu)]; aerofoil(i).y = [yl,fliplr(yu)];
aerofoil(i).r = r; aerofoil(i).yc = yc;
end
%% plot all airfoils (un-rotated)
clf; hold on
for i=1:length(aerofoil)
x = aerofoil(i).x; y = aerofoil(i).y;
r = aerofoil(i).r; yc = aerofoil(i).yc;
plot(x,y, r, yc, '--')
% plot(xu,yu,xl,yl,r,yc);
end
hold off
xlabel('x'),ylabel('y')
title( ['NACA Airfoil ', NACA]);
grid on; axis equal
%% --------------------------------------------------------------------------------------%
Yc = 0 ; % always rotate about ycenter = 0;
% Step 2) Rotate the data
angles =[1.047197551 0.756112778 0.546578176 0.413499657 0.327872784 ...
0.269927858 0.228682627 0.198042808 0.174475668 0.155828787 0.140728889 ...
0.128263734 0.117805904 0.108910765 0.101254972 0.094598036];
clf;
cmap = winter(length(angles));
hold on
for i=length(aerofoil):-1:1
% Given data vectors X and Y.
X = aerofoil(i).x; Y = aerofoil(i).y;
r = aerofoil(i).r; yc = aerofoil(i).yc;
% Want to rotate the data through angle "ang" about rotation center Xc, Yc
% Specify the coordinates of the center of rotation
Xc = 0.25*r(end) ; % Rotate about the 1/4 chord point
% The data is roated in a three-step process
% Step 1) Shift the data to the rotation center
Xs = X - Xc; % shifted data
Ys = Y - Yc;
ang = -angles(i); % in radians ??
Xsr = Xs.*cos(ang) + Ys.*sin(ang); % shifted and rotated data
Ysr = -Xs*sin(ang) + Ys*cos(ang); %
% Step 3) Un-shift the data (back to the original coordinate system)
Xr = Xsr + Xc; % Rotated data
Yr = Ysr + Yc;
plot(Xr, Yr, 'color',cmap(i,:),'linewidth',1);
h=patch(Xr, Yr, cmap(i,:),'FaceAlpha',0.2);
end % for
hold off
axis equal
grid