How to smooth an airfoil given through a set of points ?

Hi everyone!

I've written a simple code to study an airfoil geometry. The problem arises when I try to zoom in the leading edge; infact the result is a curve not smooth:

I've tryed to use interp1 but I've had some diffuculties. The code is the following one:

clc; close all; clear all;

%Import points
A = importdata ('DAE-21.txt');  %here we have 80 points starting from trailing edge, counter-clockwise direction, and returning to t.e.
x = A (:,1);  
y = A (:,2);  


do = importdata ('DAE-21_dorsoDraw.txt');
ve = importdata ('DAE-21_ventreDraw.txt');
n = length(do);

yv=interp1(ve(:,1),ve(:,2),do(:,1));      

for i=1:n-1
    X(i) = (do(i,1)+ ve(i,1))/2;
    Y(i) = (do(i,2)+ ve(i,2))/2;
end

%mean line
lm=(yv+do(:,2))/2;

figure (1)
plot (x,y,'o')
axis equal
axis ([-0.05 1.05 -0.1 0.2])
xlabel ('x/c')
ylabel ('z/c')
grid on

%Disegno del profilo: con la linea media
figure (2)
plot(x,y,'k',do(:,1),lm,'k--','LineWidth',1.2);
axis equal
axis ([-0.05 1.05 -0.1 0.2])
xlabel ('x/c')
ylabel ('z/c')
grid on

% LEADING EDGE ZOOM

figure (3)
plot (x,y,'k','LineWidth',1.2)
axis equal
axis ([-0.01 0.05 -0.03 0.05])
xlabel ('x/c')
ylabel ('z/c')
grid on

%trailing edge zoom

figure (4)
plot (x,y,'k','LineWidth',1.2)
axis equal
axis ([0.9 1.05 -0.03 0.04])
xlabel ('x/c')
ylabel ('z/c')
grid on

Whats commands can I use to smooth the airfoils or at least the leading edge ?

Answer 1

[John D'Errico]

You cannot fit a spline directly to data of that form. (I don't have your data, so I cannot give an example using it.) But for example:

theta = linspace(-pi/2,pi/2,15)';
x = cos(theta);
y = sin(theta);

what happens when we just use a spline, or interp1, for that matter on the (x,y) curve? It fails, and fails, miserably! I'll even use the method David suggested.

airfoil=fit(x,y,'smoothingspline');
plot(airfoil)

Why did it fail so badly? Becua the (x,y) pairs here are not in the form of a single valued function, that is, y(x). So for ANY x, we need to see a single value for y. But in your airfoil, you don't have that. You have a relationship that is essentially multi-valued.

So what happens is interp1 internally SORTS your data on x, then assumes y is a single valued function of x. NOT TRUE FOR YOU. And that will cause it to fail in unfortunate ways.

The common answer, IF you have a parametric form, is to use that. So, if we were to fit x(theta) and y(theta), so both functions of a third parameter that drives both x and y, then we can make it work. But you won't have that, or at least, you probably won't have it.

So what can you do? You can download my interparc function from the File Exchange.

https://www.mathworks.com/matlabcentral/fileexchange/34874-interparc

xyint = interparc(100,x,y,'spline');
plot(x,y,'bo',xyint(:,1),xyint(:,2),'r-')

If you have the curve fitting toolbox, you can use cscvn, which is able to handle a set of points that lie along a completely general path in the plane like this.

fnplt(cscvn([x,y ]'))

[Giuseppe]

Fantastic! Thank you very much, John!

Just for my curiosity: I've found another function (from file exchange)that seems to do the same things; the function is the following one:https://it.mathworks.com/matlabcentral/fileexchange/69055-interpclosed.

could you explain the differences to me ?

Best regards,

Giuseppe.

[Santiago Benito]

Hi Giuseppe,

These functions do almost the same work with a couple of exceptions.

• interpclosed works in a similar way to cscvn because it can take repeated points as inputs. If this option is used, those repeated points denote where the derivatives do not have to be equal at each side.
• interpclosed can compute the area, perimeter and centroid of the fit.

interparc does not possess these features as far as I know.

[Giuseppe]

Thank you very much, Benito! Have a good day.

[Giuseppe]

Excuse me for further question: Is it possible to increase the number of points only in a portion of my closed curve ?

p.s. I'm interested to smooth very well the zones with high curvature!

Thanks.

[Clément Carré]

Hello Giuseppe, have you found a way to increase the number of points in the portions with high curvature ?

Thanks.

Answer 2

[David Hill]

airfoil=fit(x,y,'smoothingspline');
plot(airfoil);

[John D'Errico]

David, I explain why that fails in my answer.

[David Hill]

Thanks.