MATLAB: How to smooth an airfoil given through a set of points

aerodynamicsairfoil drawingsmoothing problemtrailing edge smoothed

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 ?

Best Answer

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 ]'))

Related Solutions

MATLAB: How to move points in a curve in the normal direction

Try to reverse the sign before tol

x_n(2, 1) = x_airfoil(i, 1) + +1*tol*(y_airfoil(i+1, 1)-y_airfoil(i, 1))/sqrt(((x_airfoil(i+1, 1)-x_airfoil(i, 1))^2+(y_airfoil(i+1, 1)-y_airfoil(i, 1))^2));
y_n(2, 1) = y_airfoil(i, 1) + -1*tol*(x_airfoil(i+1, 1)-x_airfoil(i, 1))/sqrt(((x_airfoil(i+1, 1)-x_airfoil(i, 1))^2+(y_airfoil(i+1, 1)-y_airfoil(i, 1))^2));

Or the opposite

x_n(2, 1) = x_airfoil(i, 1) + -1*tol*(y_airfoil(i+1, 1)-y_airfoil(i, 1))/sqrt(((x_airfoil(i+1, 1)-x_airfoil(i, 1))^2+(y_airfoil(i+1, 1)-y_airfoil(i, 1))^2));
y_n(2, 1) = y_airfoil(i, 1) + +1*tol*(x_airfoil(i+1, 1)-x_airfoil(i, 1))/sqrt(((x_airfoil(i+1, 1)-x_airfoil(i, 1))^2+(y_airfoil(i+1, 1)-y_airfoil(i, 1))^2));

MATLAB: Extract points from spline

Ok, so you are drawing hull sections. For each z, I gather what you mean by this is you have a curve for each value of z, so some implicit function as pairs of (x,y). In your last question, z was not a factor. So I can only assume that you are trying to loft a hull from these sections.

Now, you wish to interpolate this relationship. Here is the hull section you posted from your last question.

plot(x,y,'o')
xlabel X
ylabel Y
grid on

Note that in that file you posted, we had z and x as identical vectors. But I can only read what you have said here, and try to infer what you are doing and what you really want t odo. (At least I got you out of the idea of trying to use a polynomial model to interpolate those points.)

A problem you will face in trying to interpolate these sections is there is a bit of noise in your data. As well, you have an infinite slope in this, down near x=y=1. So the third point was apparently misread, as 1.0078. My guess is it was really 1.078, not 1.0078.

[x,y]
ans =
            1            1
         1.01         2.13
       1.0078        3.261
       1.1205        4.385
       1.3478        5.492
        2.182        7.586
        2.817         8.52
         3.58        9.353
        4.429       10.099
        6.265        11.42
        7.225        12.01
        8.203        12.58
        9.188        13.13
        10.18        13.68
        12.19        14.72
        14.24        15.66
        16.36        16.44
        18.56        16.95
        20.82         17.1

Anyway, now you are asking how to interpolate that curve, trying to gain x as a function of y, for a fixed value of z? Or is your hull section actually allowing x, y, and z to vary all in one curve?

A simple way to solve the problem, based on this data, as long as the y data is monotonic (as it is here) is to use pchip. We can just interpolate at any y. So, actually, interp1 can be used, or you can use pchip, and then ppval.

yev = linspace(min(y,max(y),500);
xev = interp1(y,x,yev,'pchip');

Or do this:

yev = linspace(min(y,max(y),500);
spl = pchip(y,x);
xev = ppval(spl,yev);

Note that I carefully used y and x in that order as I called pchip, or interp1. The above will work as long as the curve at the top end is still single valued, so we have a nice functional relationship x(y). So, here we swap axes, to get the curve drawn as

plot(y,x,'o')
xlabel Y
ylabel X
grid on

So the curveas drawn here does not have an infinite slope at the ight end. And even the wiggle caused by the inaccurate third value of y will cause no problem. If you have some other desired set of points yev, then create them as you desire. So, hereif we want x at y==10, the result is as simple as:

spl = pchip(y,x);
ppval(spl,10)
ans =
       4.3074

So you can create a nice well-behaved curve, thus a list of values xev, at some fixed set of values yev.

Note that I was careful to suggest the use of pchip here, as the spline utility can cause problems if the curves get more wonky than you have shown. True splines can generate ringing-like oscillations.

As well, while I originally suggested in your last question that you could use interparc to draw the curve, you don't actually need it for this specific problem. You can get away with simply pchip, something which is already in MATLAB. One problem when you ask a question where we don't really know where you are going and what you really need to do, is it can be hard to give the best advice. In order to do that, we really need to understand your end goal.

My only fear based on your question here is that now your lofted curve actually varies in x, y, and z, all three variables moving together. If that is true, then you need to explain that to be the case. (And give me some data I can use for an example to teach you how to work with it.)

Best Answer

Related Solutions

MATLAB: How to move points in a curve in the normal direction

MATLAB: Extract points from spline

Related Question