How to calculate radiuses of an airfoil from its coordinates?

Hi,

I am trying to find the radiuses of an airfoil. I extracted the data points from a design software (Catia) to Excel. Now, l have the points of the airfoil curve, but it is defined as splines. I need to redefine the airfoil as circles.

Any help would be great.

Answer 1

[Mathieu NOE]

hello

find below a demo code that computes neutral line and curvature / radiuses of airfoil . there is also a code section that fit a circle by the Taubin method (at least 4 to 5 points are needed, whereas the curvature is computed on 3 points)

required functions are provided in attachment

hope it helps !

method with cercle fit :

Radius scatter plot the R value is color coded

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');
x = [flipud(u(:,2));flipud(l(:,2))];
y = [flipud(u(:,3));flipud(l(:,3))];
n = numel(x);% number of points

%% compute neutral axis
% centroids 
polyin = polyshape(x,y);
[xc,yc] = centroid(polyin);

% compute moment of inertia
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;

% % compute (ellipse) semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);

% axes angle
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
theta_deg = 180/pi*theta;

% % section below from link : 
% % https://leancrew.com/all-this/2018/01/python-module-for-section-properties/
% %'Principal moments of inertia (I1 I2) and orientation.'
% I_avg = (Ixx + Iyy)/2;
% I_diff = (Ixx - Iyy)/2;
% I1 = I_avg + sqrt(I_diff^2 + Ixy^2);
% I2 = I_avg - sqrt(I_diff^2 + Ixy^2);
% theta2 = atan2(-Ixy, I_diff)/2;
% theta_deg2 = 180/pi*theta2;

% draw the neutral line 
slope = tan(theta);
xl = [min(x) max(x)];
yl = yc + slope*(xl-xc);
plot(x,y,'r*-',xc,yc,'dk',xl,yl,'b--');

%% circle fit 
% let's assume we want to fit a circle to the front portion (leading edge)
ind = x<0.008;

% circle fit here then plot
par = CircleFitByTaubin([x(ind) y(ind)]);

%     Output: Par = [a b R] is the fitting circle:
%     center (a,b) and radius R
xc = par(1);
yc = par(2);
Re = par(3);

% display results in command window 
disp(['----------------------------']);
disp([' Re =  ' num2str(Re) ':']);
disp([' xc = ' num2str(xc) ':']);
disp([' yc = ' num2str(yc) ':']);
disp(['----------------------------']);


% reconstruct circle from data
n=100;
th = (0:n)/n*2*pi;
xe = Re*cos(th)+xc;
ye = Re*sin(th)+yc;

hold on 
plot(x(ind), y(ind),'db','MarkerSize',10)
plot(xe,ye,'k--')
title(' measured fitted circles')
text(xc-Re*0.5,yc + 0.75*Re,sprintf('center (%g , %g );  R=%g',xc,yc,Re))
hold off
axis equal

%% curvature computation 

[L2,R2,K2] = curvature([x y]);
% figure;
% plot(L2,R2)
% title('Curvature radius vs. cumulative curve length')
% xlabel L
% ylabel R
figure;
h = plot(x,y); grid on; axis equal
set(h,'marker','.');
xlabel x
ylabel y
title('2D curve with curvature vectors')
hold on
quiver(x,y,K2(:,1),K2(:,2));
hold off

figure
scatter(x,y,25,R2,'filled'); 
title('2D curve with radius values ')
caxis([0 3])
colormap(jet)
colorbar

[piston_pim_offset]

Thank you @Mathieu NOE ! That's close to what l need. It gives a circle for the leading edge, l'm trying to implement it for the whole airfoil.

[Mathieu NOE]

my pleasure

I assume you want to show the circles for some specific areas of the airfoil not at every sample location

[piston_pim_offset]

I want to recreate the airfoil with circles. That's what l exactly need.

[Mathieu NOE]

I guess this is what you want to do ?

here the circles do not represent the airfoil curvature

[Mathieu NOE]

then maybe this is the code you're looking for :

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');
x = [flipud(u(:,2));flipud(l(:,2))];
y = [flipud(u(:,3));flipud(l(:,3))];
nk = numel(x);% number of points

%% compute camber line
ycl = (flipud(u(:,3)) + (l(:,3)))/2; % mean of upper and lower y data
xcl = flipud(u(:,2));

plot(x, y,'.-b')
hold on 
plot(xcl, ycl,'--b')

axis equal

%% plot circles 
% first circle is obtained bit fitting 
 % let's assume we want to fit a circle to the front portion (leading edge)
ind = x<0.01*max(x);

% circle fit here then plot
par = CircleFitByTaubin([x(ind) y(ind)]);

%     Output: Par = [a b R] is the fitting circle:
%     center (a,b) and radius R
xc = par(1);
yc = par(2);
Re = par(3);

% display results in command window 
disp(['----------------------------']);
disp([' Re =  ' num2str(Re) ':']);
disp([' xc = ' num2str(xc) ':']);
disp([' yc = ' num2str(yc) ':']);
disp(['----------------------------']);


% reconstruct circle from data
nk=100;
th = (0:nk)/nk*2*pi;
xe = Re*cos(th)+xc;
ye = Re*sin(th)+yc;

plot(x(ind), y(ind),'*r','MarkerSize',10)
plot(xe,ye,'k')
title(' measured fitted circles')
% text(xc-Re*0.5,yc + 0.75*Re,sprintf('center (%g , %g );  R=%g',xc,yc,Re))

% next circles
dx = 0.05;
init = find(xcl>xc+dx,1,'first');

xcc = xcl(init):dx:xcl(end);
ycc = interp1(xcl,ycl,xcc);

for k = 1:numel(xcc)
    % find minimum radius
    Re = sqrt(min((x - xcc(k)).^2 + (y - ycc(k)).^2));
    xe = Re*cos(th)+xcc(k);
    ye = Re*sin(th)+ycc(k);

    plot(xe,ye,'k')
end
    
hold off
axis equal

[piston_pim_offset]

Not exactly, l need something like this:

[Mathieu NOE]

try this

we take for example 4 points after the leading edge and move along the top surface and do the circle fit

same logic applies for the circles at the trailing side , but for those circles we ue the bottom side data for the circle fit

in the code dx (=1) represents the jump in sample number for the next circle that is the samples used for the next circles are simply the samples of the previous circles + 1 (1 to 4 becomes 2 to 5 etc...)

if you want more spacing between circles simply increase dx (must be integer as it's used as an index))

u = readmatrix('upper.txt');
l = readmatrix('lower.txt');
xtop = flipud(u(:,2));
xbottom = flipud(l(:,2));
ytop = flipud(u(:,3));
ybottom = flipud(l(:,3));

x = [xtop;xbottom];
y = [ytop;ybottom];

nk = numel(x);% number of points

%% compute camber line
ycl = (ytop + flipud(ybottom))/2; % mean of upper and lower y data
xcl = xtop;

plot(x, y,'.-b')
hold on 
plot(xcl, ycl,'--b')

axis equal

%% plot circles 
% first circle is obtained bit fitting 
 % let's assume we want to fit a circle to the front portion (leading edge)
ind = x<0.01*max(x);

% circle fit here then plot
par = CircleFitByTaubin([x(ind) y(ind)]);

%     Output: Par = [a b R] is the fitting circle:
%     center (a,b) and radius R
xc = par(1);
yc = par(2);
Re = par(3);

% display results in command window 
disp(['----------------------------']);
disp([' Re =  ' num2str(Re) ':']);
disp([' xc = ' num2str(xc) ':']);
disp([' yc = ' num2str(yc) ':']);
disp(['----------------------------']);


% reconstruct circle from data
nk=100;
th = (0:nk)/nk*2*pi;
xe = Re*cos(th)+xc;
ye = Re*sin(th)+yc;

plot(x(ind), y(ind),'*r','MarkerSize',10)
plot(xe,ye,'k')
title(' measured fitted circles')

%% top circles close to the leading edge
dx = 1;

for k = 1:5
    
    % select top surface points only
    ind = (1:4) + k*dx;

    par = CircleFitByTaubin([xtop(ind) ytop(ind)]);

    %     Output: Par = [a b R] is the fitting circle:
    %     center (a,b) and radius R
    xc = par(1);
    yc = par(2);
    Re = par(3);
    xe = Re*cos(th)+xc;
    ye = Re*sin(th)+yc;

    plot(xe,ye,'k')
end

%% bottom circles close to the trailing edge

for k = 1:5
    
    % select top surface points only
    ind = (1:4) + k*dx;
    par = CircleFitByTaubin([xbottom(ind) ybottom(ind)]);

    %     Output: Par = [a b R] is the fitting circle:
    %     center (a,b) and radius R
    xc = par(1);
    yc = par(2);
    Re = par(3);
    xe = Re*cos(th)+xc;
    ye = Re*sin(th)+yc;

    plot(xe,ye,'k')
end

    
hold off
axis equal

[piston_pim_offset]

@Mathieu NOE thats it. I really appreciate you taking the time to help me, you are a lifesaver.

By the way, how do l get the knowledge to this by myself?

[Mathieu NOE]

my pleasure ! glad I could be of some help !

By the way, how do l get the knowledge to this by myself?

Well I don't know, my way of learning is simply to try and if I fail I look for people having successfully managed the topic. This forum and also the File exchange is a key factor in accelerating your learning curve.

[Mathieu NOE]

hello again @piston_pim_offset

just for myself and maybe for you in the future , I improved my code for this situation :

now we can see the circles are trully tangent to the top / bottom shape and the "true" camber line (in red) is slightly different from the usual approximation (average of top / bottom) in blue. the diference is more visible close to the leading edge because of the more curved shape. With a more flat airfoil , the difference would be less noticeable.

%Coordinates of Airfoil
x = [1.00000	0.94908	0.89816	0.79638	0.69491	0.59377	0.49322	0.39328	0.29389	0.19489	0.14562	0.09650	0.07203	0.04764	0.02338	0.01143	0.00000	0.01330	0.02596	0.05095	0.07573	0.10048	0.14992	0.19932	0.29821	0.39752	0.49722	0.59742	0.69804	0.79883	0.89957	0.94979	1.00000];
y = [0.00000	0.02835	0.05669	0.11138	0.15658	0.19180	0.20853	0.20678	0.18805	0.15733	0.13473	0.10764	0.09134	0.07255	0.04976	0.03287	0.00000	-0.02457	-0.02966	-0.02934	-0.02254	-0.01473	0.00237	0.02098	0.05519	0.07642	0.08566	0.07942	0.06019	0.03596	0.01324	0.00637	0.00000];

% spline interpolate the data for better results 
nk = numel(x);% number of points
N = 301; % PLS pick an ODD number to have same number of points for upper and lower sections (computation of approximate camber line)
nx = linspace(1,nk,N);
xi = spline(1:nk,x,nx);
yi = spline(1:nk,y,nx);

% % plot
% figure 
% plot(x, y,'.-b')
% hold on 
% plot(xi, yi,'.-r')
% hold off

% ok so let's move on with the new data 
x = xi;
y = yi;
clear xi yi 

%% approximate mean camber line
% this will be refined later on
[~,indle] = min(x);
[~,indte] = max(x);
    
ycl = (y(indle:-1:indte) + y(indle:end))/2; % mean of upper and lower y data
xcl = x(indle:-1:indte);

figure 
plot(x, y,'.-b')
hold on 
plot(xcl, ycl,'--b')
axis equal

%% plot circles 
 % let's assume we want to fit a circle to the front portion (leading edge)
ind = (x<0.04) & (y<0.025);

% circle fit here then plot
par = CircleFitByTaubin([x(ind)' y(ind)']);

%     Output: Par = [a b R] is the fitting circle:
%     center (a,b) and radius R
xc = par(1);
yc = par(2);
Re = par(3);

% reconstruct circle from data
nk=360;
th = (0:nk)/nk*2*pi;
xe = Re*cos(th)+xc;
ye = Re*sin(th)+yc;

plot(x(ind), y(ind),'*r','MarkerSize',5)
plot(xe,ye,'k')
title(' Fitted circles')

% now the true camber line coordinates
xcl_true = xc;
ycl_true = yc;

% next circles
dx = 0.05;
init = find(xcl>xc+dx,1,'first');

% approximate mean camber line interpolation on new vector xcc
xcc = xcl(init):dx:xcl(end);
ycc = interp1(xcl,ycl,xcc);

for k = 1:numel(xcc)-1
    % xtop, ytop
        ind = y>ycc(k); 
        xtop = x(ind);
        ytop = y(ind);    
    % xbottom ybottom
        xbot = x(~ind);
        ybot = y(~ind);
    % find nearest point of top surface vs (raw) camber line  point xcc(k),ycc(k) 
    [dtop,indtop] = min(sqrt((xtop - xcc(k)).^2 + (ytop - ycc(k)).^2));
    xtop_selected = xtop(indtop);
    ytop_selected = ytop(indtop);

    % find nearest point of bottom surface vs (raw) camber line  point xcc(k),ycc(k) 
    [dbot,indbot] = min(sqrt((xbot - xcc(k)).^2 + (ybot - ycc(k)).^2));
    xbot_selected = xbot(indbot);
    ybot_selected = ybot(indbot);
    
    
    % circle fit here then plot
    % pick 3 points on top and 3 points on bottom side
    xx = [xtop(indtop + (-1:1))';xbot(indbot + (-1:1))'];
    yy = [ytop(indtop + (-1:1))';ybot(indbot + (-1:1))'];
    
    par = CircleFitByTaubin([xx yy]);

    %     Output: Par = [a b R] is the fitting circle:
    %     center (a,b) and radius R
    xc = par(1);
    yc = par(2);
%     Re = par(3); % this is a bit too large, CircleFitByTaubin is here mostly used to find the best center of the circle but the R is overestimated
    % this is more accurate : 
    D = (sqrt((xbot_selected - xtop_selected).^2 + (ybot_selected - ytop_selected).^2));
    Re = D/2;

    xe = Re*cos(th)+xc;
    ye = Re*sin(th)+yc;
    plot(xe,ye,'k')

    % add next point to the list to create the true camber line
    xcl_true = [xcl_true xc];
    ycl_true = [ycl_true yc];
   
end

% now plot the true camber line
plot(xcl_true,ycl_true,'--r','linewidth',2)
hold off
axis equal

[piston_pim_offset]

Coordinates of my airfoil have some negative values, and this causes Newton-Taubin to not converge. How can l solve this?

[Mathieu NOE]

hmm, there must be another reason , CircleFitByTaubin works for both positive and negative data

can you share your data (and code if you have modified it)

[piston_pim_offset]

I wish I could share it, but the data is probably restricted. But your code works anyway. It gives the circles.

[piston_pim_offset]

My question is evolved to something a bit different. Now, l want to know the radiuses of the airfoil like this:

[Mathieu NOE]

that looks like a simple ellipse - what is the problem ?

you want to generate a similar plot like this one or which one of the proposed above ?

[piston_pim_offset]

Assume the ellipse as an airfoil. There will be some radiuses if we divide it to some portions. And each portion includes some number of points. What l am trying to find is the values of radiuses of the airfoil when it is divided into some portions.

[Mathieu NOE]

ok , so basically its one of my first code suggestion, simply I generated the ellipse x,y data first

now you can easily pick the last plot data and use them for this new task

n = 100;% number of points
theta = (0:n-1)'/n*2*pi;
x = 3*(cos(theta)+1);
y = sin(theta)+1;


%% circle fit 
% let's assume we want to fit a circle to the front portion (leading edge)
ind = x<0.1;

% circle fit here then plot
par = CircleFitByTaubin([x(ind) y(ind)]);

%     Output: Par = [a b R] is the fitting circle:
%     center (a,b) and radius R
xc = par(1);
yc = par(2);
Re = par(3);

% display results in command window 
disp(['----------------------------']);
disp([' Re =  ' num2str(Re) ':']);
disp([' xc = ' num2str(xc) ':']);
disp([' yc = ' num2str(yc) ':']);
disp(['----------------------------']);


% reconstruct circle from data
n=100;
th = (0:n)'/n*2*pi;
xe = Re*cos(th)+xc;
ye = Re*sin(th)+yc;

plot(x, y,'.-b')
hold on 
plot(x(ind), y(ind),'dr','MarkerSize',10)
plot(xe,ye,'k--')
title(' measured fitted circles')
hold off
axis equal

%% curvature computation 

[L2,R2,K2] = curvature([x y]);

figure;
h = plot(x,y); grid on; axis equal
set(h,'marker','.');
xlabel x
ylabel y
title('2D curve with curvature vectors')
hold on
quiver(x,y,K2(:,1),K2(:,2));
hold off
axis equal

figure
scatter(x,y,25,R2,'filled'); 
title('2D curve with radius values ')
caxis([0 3])
colormap(jet)
colorbar
axis equal