翼型轮廓几何绘图答案

【问题标题】：airfoil profile geometry plotting翼型轮廓几何绘图
【发布时间】：2014-09-21 11:50:36
【问题描述】：

我在计算坐标方面有一点问题。给出了两个列表中的翼型轮廓，具有以下示例坐标：

例子：

x_Coordinates = [1,   0.9,   0.7,   0.5,   0.3,   0.1, 0.0, ...] 
y_Coordinates = [0, -0.02, -0.06, -0.08, -0.10, -0.05, 0.0, ...]

图 1：

关于个人资料的唯一已知信息是上面的列表和以下事实：

第一个坐标始终是后沿，在上面的示例中为 (x=1, y=0)
坐标始终在前缘的底部/下侧运行，在上面的示例中位于 (0,0) 并从那里回到后缘
配置文件未标准化，可以以旋转形式存在

现在我要确定

前沿和
弧线。

直到现在，我一直使用最小的 x 坐标作为前沿。但是，这在下面是示例剖面，因为最小的 x 坐标位于剖面的上表面。

图2：

有人知道如何轻松计算/确定这些数据吗？

编辑

一个完整的样本数组数据

(1.0, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.25, 0.2, 0.15, 0.1, 0.075, 0.05, 0.025, 0.0125, 0.005, 0.0, 0.005, 0.0125, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0)
(0.00095, 0.00605, 0.01086, 0.01967, 0.02748, 0.03423, 0.03971, 0.04352, 0.04501, 0.04456, 0.04303, 0.04009, 0.03512, 0.0315, 0.02666, 0.01961, 0.0142, 0.0089, 0.0, -0.0089, -0.0142, -0.01961, -0.02666, -0.0315, -0.03512, -0.04009, -0.04303, -0.04456, -0.04501, -0.04352, -0.03971, -0.03423, -0.02748, -0.01967, -0.01086, -0.00605, -0.00095)

【问题讨论】：

如果翼型底部没有对应点，计算弧线有点棘手。你的数据来源是什么？你能假设前沿是二阶导数最大的地方吗？
Aso 是否可以获得完整的样本数组数据？
在前沿点，curvature is greatest 是真的吗？
@joojaa 谢谢帮助。我在帖子中添加了一些示例数据。我的数据来源是cpacs files。这些文件将由多种工具创建。
看起来配置文件数据实际上存储在Bernstein polynomial splines。您可以使用TiGL library 访问和评估这些曲线吗？

标签： python plot geometry curves

【解决方案1】：

好吧，我用翅膀做了很多年。

我没有任何歪斜的机翼数据，因为在您的图片上，我发现的最接近的是：

前沿不适合非平凡的翅膀

只需找到dx 的符号翻转并计算的点即可
```
dx(i)=x(i)-x(i-1)
```
然后标记dx 为正或负的区域，并找到它们之间的中间（通常dx==0 用于该区域）。将边缘点标记为ix1
弧线

对于精确的几何图形，您需要从每一侧投射法线的交点，因此：
- 从轮廓点开始i
- 在机翼内投射法线
- 搜索对面。
- 找点，法线与相反的法线相交并将两条法线划分为相同的距离
这是可行的，但非常复杂
近似弧线

不太精确的方式，但要快得多：
1. 从轮廓点开始i
2. 在对面找到离它最近的点
3. 计算它们之间的中点并将其存储为不准确的axis0 点。对所有点执行此操作i=(0-ix1)（红线）
4. 做同样的事情，但从对面的商店开始，axis1（深红色）
5. 完成后，只需找到axis0,axis1 之间的平均值
这可以通过相同的方式完成，结果是蓝色轴折线

C++ 源代码：

    List<double> pnt;   // outline 2D pnts = {x0,y0,x1,y1,x2,y2,...}
    List<double> axis;  // axis line 2D pnts = {x0,y0,x1,y1,x2,y2,...}
    int ix0,ix1;        // edge points

void compute()
    {
    int i,i0,i1;
    double d,dd;
    double *p0,*p1,*p2;
    double x0,x1,y0,y1;
    List<double> axis0,axis1;

    // find leading edge point
    ix0=0; ix1=0;
    for (p0=pnt.dat,p1=p0+2,p2=p1+2,i=2;i<pnt.num;i+=2,p0=p1,p1=p2,p2+=2)
     if ((p1[0]-p0[0])*(p2[0]-p1[0])<=0.0) { ix1=i; break; }
    // find axis0: midpoint of i0=(0-ix1) i1=find closest from (ix1,pnt.num)
    for (i0=2,i1=pnt.num-2;i0<ix1-2;i0+=2)
        {
        x0=pnt[i0+0];
        y0=pnt[i0+1];
        for (d=-1.0,i=i1;i>ix1+2;i-=2)
            {
            x1=pnt[i+0];
            y1=pnt[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||(dd<=d)) { i1=i; d=dd; }
            }
        if (d>=0.0)
            {
            x1=pnt[i1+0];
            y1=pnt[i1+1];
            axis0.add(0.5*(x0+x1));
            axis0.add(0.5*(y0+y1));
            }
        }
    // find axis1: midpoint of i0=(ix1,pnt.num) i1=find closest from (0,ix1)
    for (i1=2,i0=pnt.num-2;i0>ix1+2;i0-=2)
        {
        x0=pnt[i0+0];
        y0=pnt[i0+1];
        for (d=-1.0,i=i1;i<ix1-2;i+=2)
            {
            x1=pnt[i+0];
            y1=pnt[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||(dd<=d)) { i1=i; d=dd; }
            }
        if (d>=0.0)
            {
            x1=pnt[i1+0];
            y1=pnt[i1+1];
            axis1.add(0.5*(x0+x1));
            axis1.add(0.5*(y0+y1));
            }
        }
    // find axis: midpoint of i0=<0-axis0.num) i1=find closest from <0-axis1.num)
    axis.add(pnt[ix0+0]);
    axis.add(pnt[ix0+1]);
    for (i0=0,i1=0;i0<axis0.num;i0+=2)
        {
        x0=axis0[i0+0];
        y0=axis0[i0+1];
        for (d=-1.0,i=i1;i<axis1.num;i+=2)
            {
            x1=axis1[i+0];
            y1=axis1[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||(dd<=d)) { i1=i; d=dd; }
            }
        if (d>=0.0)
            {
            x1=axis1[i1+0];
            y1=axis1[i1+1];
            axis.add(0.5*(x0+x1));
            axis.add(0.5*(y0+y1));
            }
        }
    axis.add(pnt[ix1+0]);
    axis.add(pnt[ix1+1]);
    }

List<double> xxx; 只是我的动态列表模板，和double xxx[]; 一样
xxx.add(5);将 5 添加到列表末尾
xxx[7]访问数组元素
xxx.num 是数组的实际使用大小
xxx.reset() 清空数组并设置 xxx.num=0

[edit1] 正确的前沿点

有一个疯狂的想法来找到运行中的边缘点加上一些代码调整，结果对我来说已经足够好了:) 所以首先解释一下：

轴的算法保持不变，但不是ix1 bound 只使用尚未使用的点...如果没有找到停止点（上图情况），也只计算有效的最近点（在对面）。从这一点找到距离最后一个轴点最远的点，这是前缘点。

这种方法的输出非常准确（axis0,axis1 更接近）

现在是 C++ 代码：

void compute()
    {
    int i,i0,i1,ii,n=4;
    double d,dd;
    double x0,x1,y0,y1;
    List<double> axis0,axis1;
    ix0=0; ix1=0;

    // find axis0: midpoint of i0=(0-ix1) i1=find closest from (ix1,pnt.num)
    for (i0=0,i1=pnt.num-2;i0+n<i1;i0+=2)
        {
        x0=pnt[i0+0];
        y0=pnt[i0+1];
        i=i1+n; if (i>pnt.num-2) i=pnt.num-2; ii=i1;
        for (d=-1.0;i>i0+n;i-=2)
            {
            x1=pnt[i+0];
            y1=pnt[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||((dd<=d)&&(dd>1e-10))) { i1=i; d=dd; }
            if ((d>=0.0)&&(dd>d)) break;
            }
        if (d>=0.0)
            {
            if (i1-i0<=n+2) { i1=ii; break; } // stop if non valid closest point found
            x1=pnt[i1+0];
            y1=pnt[i1+1];
            axis0.add(0.5*(x0+x1));
            axis0.add(0.5*(y0+y1));
            }
        }
    // find leading edge point (the farest point from last found axis point)
    x0=axis0[axis0.num-2];
    y0=axis0[axis0.num-1];
    for (d=0.0,i=i0;i<=i1;i+=2)
        {
        x1=pnt[i+0];
        y1=pnt[i+1];
        dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
        if (dd>d) { ix1=i; d=dd; }
        }
    axis0.add(pnt[ix1+0]);
    axis0.add(pnt[ix1+1]);

    // find axis1: midpoint of i0=(ix1,pnt.num) i1=find closest from (0,ix1)
    for (i1=0,i0=pnt.num-2;i0+n>i1;i0-=2)
        {
        x0=pnt[i0+0];
        y0=pnt[i0+1];
        i=i1-n; if (i<0) i=0; ii=i1;
        for (d=-1.0;i<i0-n;i+=2)
            {
            x1=pnt[i+0];
            y1=pnt[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||((dd<=d)&&(dd>1e-10))) { i1=i; d=dd; }
            if ((d>=0.0)&&(dd>d)) break;
            }
        if (d>=0.0)
            {
            if (i0-i1<=n+2) { i1=ii; break; } // stop if non valid closest point found
            x1=pnt[i1+0];
            y1=pnt[i1+1];
            axis1.add(0.5*(x0+x1));
            axis1.add(0.5*(y0+y1));
            }
        }
    // find leading edge point (the farest point from last found axis point)
    x0=axis1[axis1.num-2];
    y0=axis1[axis1.num-1];
    for (d=0.0,i=i1;i<=i0;i+=2)
        {
        x1=pnt[i+0];
        y1=pnt[i+1];
        dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
        if (dd>d) { ix1=i; d=dd; }
        }
    axis1.add(pnt[ix1+0]);
    axis1.add(pnt[ix1+1]);

    // find axis: midpoint of i0=<0-axis0.num) i1=find closest from <0-axis1.num)
    for (i0=0,i1=0;i0<axis0.num;i0+=2)
        {
        x0=axis0[i0+0];
        y0=axis0[i0+1];
        for (d=-1.0,i=i1;i<axis1.num;i+=2)
            {
            x1=axis1[i+0];
            y1=axis1[i+1];
            dd=((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0));
            if ((d<0.0)||(dd<=d)) { i1=i; d=dd; }
            }
        if (d>=0.0)
            {
            x1=axis1[i1+0];
            y1=axis1[i1+1];
            axis.add(0.5*(x0+x1));
            axis.add(0.5*(y0+y1));
            }
        }
    }

常量n=4 只是为了安全重叠搜索最近的点，它应该是pnt.num 的一小部分。有时最近点在最后找到的最近点之前，这取决于展位侧面的曲率。太大的n 会导致速度变慢，如果n>pnt.num/4 也可能导致输出无效。

如果太小，那么对于较小的曲率半径会降低精度，这种方法取决于足够的点覆盖。如果机翼的采样点数太少，可能会导致不准确。源代码几乎是您可以选择要记住的 ix1 的 3 倍（从第一次或第二次搜索中）它们是相邻点

测试配置文件：

1.000000 0.000000
0.990000 0.006719
0.980000 0.013307
0.970000 0.019757
0.960000 0.026064
0.950000 0.032223
0.940000 0.038228
0.930000 0.044075
0.920000 0.049759
0.910000 0.055276
0.900000 0.060623
0.890000 0.065795
0.880000 0.070790
0.870000 0.075604
0.860000 0.080234
0.850000 0.084678
0.840000 0.088935
0.830000 0.093001
0.820000 0.096876
0.810000 0.100558
0.800000 0.104046
0.790000 0.107339
0.780000 0.110438
0.770000 0.113342
0.760000 0.116051
0.750000 0.118566
0.740000 0.120887
0.730000 0.123016
0.720000 0.124954
0.710000 0.126702
0.700000 0.128262
0.690000 0.129637
0.680000 0.130829
0.670000 0.131839
0.660000 0.132672
0.650000 0.133331
0.640000 0.133818
0.630000 0.134137
0.620000 0.134292
0.610000 0.134287
0.600000 0.134127
0.590000 0.133815
0.580000 0.133356
0.570000 0.132755
0.560000 0.132016
0.550000 0.131146
0.540000 0.130148
0.530000 0.129030
0.520000 0.127795
0.510000 0.126450
0.500000 0.125000
0.490000 0.123452
0.480000 0.121811
0.470000 0.120083
0.460000 0.118275
0.450000 0.116392
0.440000 0.114441
0.430000 0.112429
0.420000 0.110361
0.410000 0.108244
0.400000 0.106085
0.390000 0.103889
0.380000 0.101663
0.370000 0.099414
0.360000 0.097148
0.350000 0.094870
0.340000 0.092589
0.330000 0.090309
0.320000 0.088037
0.310000 0.085779
0.300000 0.083541
0.290000 0.081329
0.280000 0.079149
0.270000 0.077006
0.260000 0.074906
0.250000 0.072855
0.240000 0.070858
0.230000 0.068920
0.220000 0.067047
0.210000 0.065242
0.113262 0.047023
0.110002 0.042718
0.106385 0.038580
0.102428 0.034615
0.098146 0.030832
0.093556 0.027239
0.088673 0.023844
0.083516 0.020652
0.078101 0.017670
0.072448 0.014904
0.066574 0.012361
0.060499 0.010044
0.054241 0.007958
0.047820 0.006108
0.041256 0.004497
0.034569 0.003129
0.027779 0.002005
0.020907 0.001129
0.013972 0.000502
0.006997 0.000126
0.000000 0.000000
0.000000 0.000000
-0.003997 0.000126
-0.007972 0.000502
-0.011907 0.001129
-0.015779 0.002005
-0.019569 0.003129
-0.023256 0.004497
-0.026820 0.006108
-0.030241 0.007958
-0.033499 0.010044
-0.036574 0.012361
-0.039448 0.014904
-0.042101 0.017670
-0.044516 0.020652
-0.046673 0.023844
-0.048556 0.027239
-0.050146 0.030832
-0.051428 0.034615
-0.052385 0.038580
-0.053002 0.042718
-0.053262 0.047023
-0.053153 0.051484
-0.052659 0.056093
-0.051768 0.060841
-0.050467 0.065717
-0.048744 0.070711
-0.046588 0.075813
-0.043988 0.081012
-0.040935 0.086297
-0.037420 0.091658
-0.033435 0.097082
-0.028972 0.102558
-0.024025 0.108074
-0.018589 0.113618
-0.012657 0.119178
-0.006228 0.124741
0.000704 0.130295
0.008139 0.135828
0.016079 0.141326
0.024525 0.146777
0.033475 0.152169
0.042930 0.157488
0.052885 0.162722
0.063339 0.167858
0.074287 0.172883
0.085723 0.177784
0.097643 0.182549
0.110038 0.187166
0.122902 0.191621
0.136226 0.195903
0.150000 0.200000
0.164214 0.203899
0.178856 0.207590
0.193914 0.211059
0.209376 0.214297
0.225227 0.217291
0.241453 0.220032
0.258039 0.222509
0.274968 0.224711
0.292223 0.226629
0.309787 0.228254
0.327641 0.229575
0.345766 0.230585
0.364142 0.231274
0.382749 0.231636
0.401566 0.231662
0.420570 0.231345
0.439740 0.230679
0.459054 0.229657
0.478486 0.228274
0.498015 0.226525
0.517615 0.224404
0.537262 0.221908
0.556930 0.219032
0.576595 0.215775
0.596231 0.212132
0.615811 0.208102
0.635310 0.203684
0.654700 0.198876
0.673956 0.193679
0.693050 0.188091
0.711955 0.182115
0.730644 0.175751
0.749091 0.169002
0.767268 0.161869
0.785149 0.154357
0.802706 0.146468
0.819913 0.138207
0.836742 0.129580
0.853169 0.120591
0.869166 0.111246
0.884707 0.101553
0.899768 0.091518
0.914322 0.081149
0.928345 0.070455
0.941813 0.059445
0.954701 0.048128
0.966987 0.036514
0.978646 0.024614
0.989658 0.012439
1.000000 0.000000

【讨论】：

感谢您的解决方案。你的近似弧线的想法很棒。我目前开发了一个从 cpacs 文件中读取配置文件数据的工具。这些配置文件可以是非常规的，也可以以旋转方式存储。我的中心问题是计算此类配置文件的领先优势。我在wikipedia 上找到了重要的侧翼。因此，我目前没有此个人资料的数据。
@René 我目前正在寻找非平凡机翼的前沿点（我发现该轮廓的只有图像没有 dat 文件，它是涡轮叶片而不是机翼）已将我的旧机翼类代码移植到较新的编译器为此，我用一些 sin、cos 函数来近似翼型，现在我检查一种有趣的方法来在轴搜索过程中找到边缘点。有趣的是，与我现在的代码风格相比，10 多年的代码看起来如何 :)
@René 我做到了...在我的答案中添加了edit1。我知道这有点粗略，但它也可以添加我的近似配置文件
哇，这太酷了，一开始它似乎效果很好:-) 非常感谢。
@René 它适用于几乎所有我的个人资料。我发现只有一个不正常（一个搜索正常，第二个不正常），但这几乎可以肯定是配置文件中的一个错误，因为与它非常相似的配置文件可以正常工作（我的旧配置文件在旧的飞机自动测量系统中非常有问题）数据很不稳定，需要一些平滑（这就是我的原始代码，顺便说一句，这是我的第一个商业产品：））......很高兴能帮上忙......

【解决方案2】：

有点过时了，但我还是看到了这篇文章。

我对前缘圆的解决方案是首先通过翼型坐标插入一条样条线。这给出了翼型的平滑参数表示。然后计算参数曲线的曲率和曲率圆，取最小的圆。这定义了前缘位置并返回那里的半径。

下面有两个 Python 函数（依赖于 numpy 和 scipy 包）这样做：

def spline(self, x, y, points=200, degree=2, evaluate=False):
    """Interpolate spline through given points

    Args:
        spline (int, optional): Number of points on the spline
        degree (int, optional): Degree of the spline
        evaluate (bool, optional): If True, evaluate spline just at
                                   the coordinates of the knots
    """

    # interpolate B-spline through data points
    # returns knots of control polygon
    # tck ... tuple (t,c,k) containing the vector of knots,
    # the B-spline coefficients, and the degree of the spline.
    # u ... array of the parameters for each knot
    # NOTE: s=0.0 is important as no smoothing should be done on the spline
    # after interpolating it
    tck, u = interpolate.splprep([x, y], s=0.0, k=degree)

    # number of points on interpolated B-spline (parameter t)
    t = np.linspace(0.0, 1.0, points)

    # if True, evaluate spline just at the coordinates of the knots
    if evaluate:
        t = u

    # evaluate B-spline at given parameters
    # der=0: returns point coordinates
    coo = interpolate.splev(t, tck, der=0)

    # evaluate 1st derivative at given parameters
    der1 = interpolate.splev(t, tck, der=1)

    # evaluate 2nd derivative at given parameters
    der2 = interpolate.splev(t, tck, der=2)

    spline_data = [coo, u, t, der1, der2, tck]

    return spline_data

计算参数曲线曲率属性的函数：

    def getCurvature(spline_data):
    """Curvature and radius of curvature of a parametric curve

    der1 is dx/dt and dy/dt at each point
    der2 is d2x/dt2 and d2y/dt2 at each point

    Returns:
        float: Tuple of numpy arrays carrying gradient of the curve,
               the curvature, radiusses of curvature circles and
               curvature circle centers for each point of the curve
    """

    coo = spline_data[0]
    der1 = spline_data[3]
    der2 = spline_data[4]

    xd = der1[0]
    yd = der1[1]
    x2d = der2[0]
    y2d = der2[1]
    n = xd**2 + yd**2
    d = xd*y2d - yd*x2d

    # gradient dy/dx = dy/du / dx/du
    gradient = der1[1] / der1[0]

    # radius of curvature
    R = n**(3./2.) / abs(d)

    # curvature
    C = d / n**(3./2.)

    # coordinates of curvature-circle center points
    xc = coo[0] - R * yd / np.sqrt(n)
    yc = coo[1] + R * xd / np.sqrt(n)

    return [gradient, C, R, xc, yc]

示例：

详情：

请注意，我还在前缘使用了样条线的细化。相应的算法在此不作介绍（作为题外话）。

【讨论】：