python电化学cv曲线怎么画_我如何使用参数自动绘制翼型

引言

在用遗传算法优化翼型之前，我们需要建立一个足够大的翼型种群库。这时候就需要我们参数化绘制翼型，自动生成第一批种子。

翼型的构成

首先，我们得知道翼型是怎么构成的。常用翼型测算工具的都知道，单弯度翼型有四个重要的参数，那就是相对厚度及其位置，最大弯度及其位置。

比如用xfoil导入一个翼型后，输入GDES，TSET命令就可以看到一个翼型被分解成两条曲线后的状态。（厚度曲线是对称的）

这四个参数能决定一个翼型的大致样貌，但是不能完全决定一个翼型。也就是用不同的方法生成一个翼型，即使参数相同，翼型也会有所不同。而形状上一点点的不同，会对气动有着蝴蝶效应般的影响。

如下两个翼型，一个是大名鼎鼎的MH 114，一个是上图中的foil 002。可以看出，参数（几乎）是一致的，但是形状上却有一定区别。计算上，foil 002性能远不如MH 114。

但这些都是后话，重点是如何用不同的方法把给定参数的翼型画出来。

二次曲线法理论计算

思路就是先确定弯度线和厚度线，然后两者一结合就可以画出翼型了。 我们可以再看一下第一幅图，看看这两条曲线有什么性质。可以很直观地得到，厚度曲线

经过

点，并在这一点取得极大值。其中

是最大厚度位置，

是最大相对厚度。弯度曲线同理。另外，厚度曲线在

处的切线垂直于x轴

首先，我们可以想到用一个二次的曲线去模拟翼型曲线。二次曲线的公式是

。由于翼型过

和

点，很容易知道

和

。

对

求导得：

把在

导数为零代入，即令

得（用

代替），即

。再把这一点代入曲线表达式（这一点也在曲线上），化简得：

也就是说到目前为止只有

和

是自由参数，其他都由相对厚度及其位置所决定

现在考虑在

处的导数。

对

求导得：

把

处

代入得：

。这样就只有

是个自由参数了。但是，这个参数也有一定的限制，对于函数

，我们只能保证

是其一个解，但是

时

却不一定为

，也就是会有两个解。

时同理。我们要做的就是要保证

和

是同侧解，也就是在

的取值上均取一侧值，否则图像可能变成这样：

蓝色标记的为同侧解 。

分别求出

为

或

时的另一个解，是

和

。由两者同号，得出

。而

的值与

有关。这样前期理论就完成的差不多了。代码实现

先定义一个行向量x（数组），从零到一，均分100份。输入四个参数，两两组合可以求出一条曲线（剩下的参数根据约束随意取一个）。比如厚度曲线：

c = -4

e = 0

a = c*v**2/p**2

b = a*(1-2*p)/v

d = -a

之后对每个x，也就是行向量x中的项求对应y值。但是不必写成循环，直接用广播即可。

A = c

B = b*x+e

C = a*x**2+d*x

Delta = B**2-4*A*C

y = (-B-np.sqrt(Delta))/(2*A)

前面各项就是一元二次方程的系数（x也看做常数），Delta是判别式。如果参数都是按照约束条件取的，判别式就不会为负

对于弯度曲线也是一样的，就是e的取值变得自由了。另外，根据之前的公式，曲线两端的导数均可求，左边是-a/e，右边是a/(b+e)，可以看出e的实际意义。 这样一个翼型就由四个主要参数和三个次要参数（两条曲线的c和弯度曲线的e）决定。

不过这样画出来的曲线是有问题的，因为x是等分的，所以在前缘附近由于斜率较大点会比较稀疏。所以可以根据初生成的翼型形状重新确定一下x的分布，再画一次翼型。然后将弯度曲线分别加上或减去厚度曲线，就得到最终的翼型，如图所示：

贝塞尔曲线法理论介绍

接下来介绍一下另一种方法，叫做Bezier曲线法。

这种方法的思路在于把曲线分成最大厚度或者最大弯度所在位置的前后两个部分，在这一点上满足一定的连续条件，比如切线和曲率相等。

贝塞尔曲线（CAD中输入spl->m->cv即可调用）是当今计算机辅助制造中最常用的造型曲线之一。具体的原理我就不多做介绍，重点讲一讲他的几个性质。

先给出三个控制点的曲线方程，正好也是二次的：

其中，

是第

个控制点。贝塞尔曲线由构造点确定，其中始末两个构造点在曲线上，其余不在。

两条曲线端点相连，相切连续的充要条件是

，而曲率连续的必要条件是

，充分条件还得带上相切连续。代码实现

这样，在给定

处的斜率、极值和极值点后，我们就可以把五个端点写出来了。

P = np.array([(0, 0), (val*k0, val), (pos, val), (0, 0), (1, 0)])

P[3] = (-P[0]+P[4]+2*P[1])/2

第四个点P[3]由其余四个点计算得来。

再定义一个列向量t，从零到一，用广播的思想就可以把整个曲线画出来。

t = np.linspace(0, 1, 101).reshape(101, 1)

cv0 = (1-t)**2*P[0]+2*t*(1-t)*P[1]+t**2*P[2]

cv1 = (1-t)**2*P[2]+2*t*(1-t)*P[3]+t**2*P[4]

cv = np.vstack((cv0, cv1))

需要注意的是，这里的t和之前的x不一样，不是横坐标，而是要把横坐标和纵坐标一起求出来 ，这样就造成了一个问题，弯度曲线和厚度曲线的对应点找不到了，所以不能直接相加，因此还需要对曲线进行线性插值。

for xi in x:

left = np.argmax(cv[0][cv[0]<=xi])

right = np.argmin(cv[0][cv[0]>=xi]

dista = cv[1][right]-cv[1][left]

new_cv.append((xi, cv[1][left]if dista==0 else(cv[0][right]-xi)*cv[1][left]/dista+(xi-cv[0][left])*cv[1][right]/dista))

这样翼型就画好了（这里没处理前缘点疏密问题，处理方法同上，不再赘述）：

小结

今天主要介绍了如何参数绘制翼型。共介绍了两种方法，分别是二次曲线法和贝塞尔曲线法。但是绘制翼型远不止这么两种方法，但是思路是一致的。

都看到这了，留个赞呗~