13
$\begingroup$

I am trying to produce a smooth spline curve through a set of airfoil coordinates. An example for the data I am trying to use is here, and I am reading the file like so to get rid of some extraneous stuff:

airfoildata = Drop[Take[ Drop[ Import["n64015.dat"], 3], 53], {27}]; airfoildata = Join[Reverse[Drop[Take[airfoildata, 26], 1]], Drop[airfoildata, 26]];

ListLinePlot gives a plot of the airfoil that's alright, but you can see some "kinks" in the plot near the leading edge, where data are sparse relative to the curvature. So, I tried getting a real smooth curve via B-Splines, using this code:

pts = Take[airfoildata, 20]; n = Length[pts]; dist = Accumulate[ Table[EuclideanDistance[pts[[i]], pts[[i + 1]]], {i, Length[pts] - 1}]]; param = N[Prepend[dist/Last[dist], 0]]; deg = 3; knots = Join[ConstantArray[0, deg], Range[0, 1, 1/(n - deg)], ConstantArray[1, deg]]; m = Table[BSplineBasis[{deg, knots}, j - 1, param[[i]]], {i, n}, {j, n}]; ctrlpts = LinearSolve[m, pts];

and then plotting it via

ListPlot[airfoildata, Prolog -> BSplineCurve[ctrlpts], PlotStyle -> Directive[Red, PointSize[Medium]]]

Using only the first 20 points as in the above code fragment, things look good, but if I use 21 points I can see numerical instability creeping in, for 22 points things get truly awful, and beyond that the linear system for the control points has no solution.

So here's my question: How can I obtain a smooth curve (with a nicely rounded leading edge near the origin) that respects the given airfoil coordinates? It looks like a standard B-spline will not work, and somehow this seems to be related to the fact that I have a relatively large number of points. However, I have used various CAD packages in the past that had no problems creating a smooth curve through data like the ones I have.

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ Since the data you linked corresponds to a NACA airfoil, you might be interested in an implementation of the formulae describing the profile that came up previously on the site. See: Proper use of a formula to produce a cambered airfoil. $\endgroup$ Commented Jul 27, 2016 at 20:53

2 Answers 2

Reset to default
9
$\begingroup$

After permuting the points so that the cusp is the first point (just like in Gypaets's answer), the methods of this answer can be used:

airfoil = Flatten[MapAt[Most @* Reverse, DeleteCases[Split[Drop[ Import["http://m-selig.ae.illinois.edu/ads/coord/n64015.dat"], 3], # =!= {} && #2 =!= {} &], {{}}], 1], 1]; (* Lee's method, http://dx.doi.org/10.1016/0010-4485(89)90003-1 *) parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ] tvals = parametrizeCurve[airfoil]; m = 3; (* degree of the B-spline *) (* knots for interpolating B-spline *) knots = Join[ConstantArray[0, m + 1], MovingAverage[ArrayPad[tvals, -1], m], ConstantArray[1, m + 1]]; (* basis function matrix *) bas = Table[BSplineBasis[{m, knots}, j - 1, tvals[[i]]], {i, Length[airfoil]}, {j, Length[airfoil]}]; ctrlpts = LinearSolve[bas, airfoil]; Graphics[{{Directive[AbsoluteThickness[3], ColorData[97, 1]], BSplineCurve[ctrlpts, SplineDegree -> 3, SplineKnots -> knots]}, {AbsolutePointSize[5], Point[airfoil]}}, Axes -> None, Frame -> True]

airfoil

Replace BSplineCurve[] with BSplineFunction[] if needed.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Very nice solution, answering my (implied) question on a "good" parameterization of the curve for the purpose of interpolation. Thanks! I'll mark this as my best answer. $\endgroup$ Commented Jul 29, 2016 at 17:02
  • $\begingroup$ The code in your question used chord-length parametrization (corresponding to setting a to 1 in parametrizeCurve[]), which should also work; the snag you hit into was that you did not permute the points so that the cusp is at the start and end. But Lee shows that centripetal parametrization is often a more "natural" choice for many configurations. $\endgroup$ Commented Jul 29, 2016 at 17:06
5
$\begingroup$

In an airfoil you probably want a sharp trailing edge and a rounded leading edge. Therefore I would choose the trailing edge as first and last point for your spline:

path = "http://m-selig.ae.illinois.edu/ads/coord/n64015.dat"; airfoilData = Drop[Take[Drop[Import@path,3], 53], {27}]; upperSide = airfoilData[[1 ;; Length@airfoilData/2]]; lowerSide = airfoilData[[Length@airfoilData/2 + 1 ;; -1]]; airfoilDataOrdered = Reverse@upperSide~Join~lowerSide[[2 ;; -1]];

Create spline:

f = BSplineFunction[airfoilDataOrdered];

Plot:

ListPlot[{airfoilDataOrdered, f /@ Range[0, 1, 0.01]}]

EDIT 1

To go through the points you can do a parametric interpolation.

f2 = Interpolation[Table[{i, airfoilDataOrdered[[i]]}, {i, Length@airfoilDataOrdered}], Method -> "Spline", InterpolationOrder -> 5]; ParametricPlot[f2[x], {x, 1, Length@airfoilDataOrdered}, PlotRange -> {{0, 0.03}, {-0.03, 0.03}}]

The nose should now be ok.

EDIT 2

Polar approach:

iPoints = Table[{π - ArcTan[-airfoilDataOrdered[[i, 1]] + 1/4, airfoilDataOrdered[[i, 2]]], airfoilDataOrdered[[i]] - {1/4, 0}}, {i, Length@airfoilDataOrdered}]; iPoints[[-1, 1]] = 2 π; f3 = Interpolation[iPoints, Method -> "Spline", InterpolationOrder -> 2]; ParametricPlot[f3[x], {x, 0, 2 π}, PlotRange -> All]
Improve this answer
$\endgroup$
5
  • $\begingroup$ Hmm, this kind of works, except for the fact that this takes the airfoil coordinates as the control points of the spline. As a consequence, the resulting function does not, in general, go through the given points. For example, in your solution above you can see that the leading edge (the very front of the airfoil) now sits around x=0.0025, not at x=0 as it should. $\endgroup$ Commented Jul 27, 2016 at 20:42
  • $\begingroup$ @Pirx Answer edited with a parametric interpolation. $\endgroup$ Commented Jul 27, 2016 at 22:07
  • $\begingroup$ Awesome, that's exactly what I was looking for. Thanks! $\endgroup$ Commented Jul 27, 2016 at 22:15
  • $\begingroup$ One more comment: Your solution above uses the index of the airfoil coordinates as the independent variable in a parametric definition of the airfoil coordinates. I had this idea that it might be better to preserve the information on the spacing of the data points, doing something like f3 = Interpolation[ Table[{Sign[26 - i] airfoildata[[i, 1]], airfoildata[[i]]}, {i, Length@airfoildata}], Method -> "Spline", InterpolationOrder -> 3]; It turns out this works also, but if you increase the interpolation order this approach breaks down completely. $\endgroup$ Commented Jul 27, 2016 at 22:39
  • $\begingroup$ @Pirx In that case keep it low. You will almost always have problems with higher Interpolation orders. I have also added a polar approach where the nose seems to be a bit more rounded than with your spacing. $\endgroup$ Commented Jul 28, 2016 at 7:30

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.