B-spline fit of airfoil coordinates

16 ビュー (過去 30 日間)
Marco Zanichelli
Marco Zanichelli 2021 年 9 月 3 日
Hello everybody,
I am an aerospace Engineer, I am currently working in an airfoil optimization project. I want to create an airfoil parametrization starting from a RAE2822 of which I know the coordinates.
At the beginning I was trying to get a unique B-spline representation for the whole airfoil curve, i.e. to find a unique curve that fits the airfoil curve. I have not been able to do so (does anybody have an idea?) . Anyway, I found the spap2 function and I've been using it to find the B-spline best-fit to the coordinates I have, keeping separate the lower and upper part of the airfoil. Trying some different knot sequences and found some that seem promising, but I am getting a weird behaviour at the Leading Edge (i.e. the "nose" of the airfoil). The B-spline curves follow really accurately the original curve but show a "pointy" behaviour. I am now attaching my script.
I am guessing this problem comes from the large derivative I have at the beginning of the curve, but I am not really sure about that. Can somebody help ?
Thank you in advance
filename = 'RAE2822.txt'; % attached here
coords = readtable(filename);
x_0 = coords.Var1;
y_0 = coords.Var2;
% separation of upper and lower parts of the airfoil
xaxis = sort(x_0(1:65));
x_1 = x_0(1:65);
y_1 = y_0(1:65);
x_2 = [0; x_0(66:end)];
y_2 = y_0(65:end);
% plot of the original coordinates
plot(x_1, y_1, x_2, y_2); hold on;
% definition of the knot sequence
x = [0,0.001 0.005, 0.01, 0.025, 0.05, 0.2, 0.55, 1];
degree = 3 %b - spline degree
degree = 3
knots = augknt(x,degree+1)
knots = 1×15
0 0 0 0 0.0010 0.0050 0.0100 0.0250 0.0500 0.2000 0.5500 1.0000 1.0000 1.0000 1.0000
w1 = ones(size(x_1)); %
w1(end) = 1000;
w2 = ones(size(x_2));
w2(1) = 1000;
%least squares b-spline fit of the upper curve (suction side (AKA SS))
BS_SS = spap2(knots, degree+1, x_1, y_1,w1);
%least squares b-spline fit of the lower curve (pressure side (AKA PS))
BS_PS = spap2(knots, degree+1, x_2, y_2, w2);
fnplt(BS_PS); axis equal
xlim([-0.01, 0.05]) % comment if you want to see the whole parametrization. Here you can see the problem

回答 (1 件)

Christopher McCausland
Christopher McCausland 2021 年 9 月 3 日
Hi Marco,
Thats a really cool problem!
Looking at it I would assume that the resoloution between x = 0 and x = 0.01 is not great enough, therefore an error is introduced. Looking at spap2 it is the least-squares spline approximation method so it is an approximation. Increasing the number of knots will improve the approximation at the cost of producing a more complex curve.
Consider adding more knot values between x = 0 and 0.01 to reduce this error and provide a better result.
Let me know if this helps,
  5 件のコメント
Christopher McCausland
Christopher McCausland 2021 年 9 月 6 日
Hi Marco,
Yep, more knots will improve the approximation. However I see that you want to keep it around 7 for a simplier curve which is a valid point.
In order to do so we would have to think about where we are placing the knots. Remeber that the curve is forced through the knot point and the curve between knots is approximate. Therefore if you move more of your 7 knots (3/4) to where the curve has the highest differential (x = 0 => x =~ 0.015) - rather than spaced equally along the x axis - you will get a better approximation for the tighter geometry at the cost of a worse approximation along the rest of the curve, however as the rest of the curve is flatter this worse approximation should be far less noticable.
So in short, yes it is possiable to get a better approximation with a minimum number of knots by opting to place more of the knots at the more compex geometry.
Thank you for unsupressing 'knots' too, thats what I was looking for.
I think what @darova is trying to alude to is that what you are doing is commonly refered to as cubic spline interpolation. See MATLAB help Spline, that being a said as you already have the script set up and running I wouldn't bother changing to spline. Having used the spine function before for active high pass filtering and its great!


Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by