张 骥， 朱春钢， 冯仁忠， 刘明明， 张恒洋

(1. 大连理工大学数学科学学院，辽宁 大连 116024；2. 北京航空航天大学数学与系统科学学院，北京 100083)

一种改进的B样条翼型参数化方法

张 骥1， 朱春钢1， 冯仁忠2， 刘明明1， 张恒洋1

(1. 大连理工大学数学科学学院，辽宁 大连 116024；2. 北京航空航天大学数学与系统科学学院，北京 100083)

翼型设计是空气动力学研究的一项重要内容，翼型的参数化结果将影响翼型的优化设计。为了减少翼型优化中的设计变量，保证优化结果的光滑性与C2连续，在优化过程中控制翼型几何特性的变化范围，提出了一种改进的B样条参数化方法。用一条三次非均匀B样条曲线表示翼型，翼型数据的参数化过程中主要运用了B样条曲线拟合算法，并且在一般的B样条曲线拟合算法的基础上加入了对曲线的法向约束，通过迭代得到最终的参数化结果。实验结果表明，该方法可以很好的拟合典型的翼型数据，得到的翼型参数化结果不仅光滑，满足 C2条件，而且所得翼型函数的参数个数比传统的参数化方法有了进一步的减少，更有利于之后翼型的优化设计。

翼型；B样条曲线；翼型参数化；翼型优化

翼型设计是空气动力学研究的一项重要内容，翼型作为飞行器翼面部件的截面形状，其几何外形设计对整个翼面部件甚至是整个飞行器的性能均存在重要影响。随着现代飞行器的设计要求越来越多，性能目标不断提高，在设计过程中必须对飞行器的几何外形进行精细化设计，因此对翼型的设计优化提出了更高的要求。为了缩短研发周期，降低研发费用，将基于气动性能数值计算的翼型优化设计方法与风洞实验相结合已成为目前翼型研究的发展方向[1-2]。

为了实现基于数值计算的翼型优化，需要给出相应的翼型参数化方法，即用含参数的翼型函数拟合由离散数据点表示的待优化翼型，再选择翼型函数的适当参数作为设计变量，结合优化算法与流场计算实现优化设计。在气动优化过程中，翼型参数化方法的选择除了影响优化算法类型的选择外，还影响计算时间和资源，设计空间的范围及设计空间中翼型几何外形是否光滑，设计空间中是否包含有意义的优化结果，在保持结果有足够的自由度和光滑性的前提下尽可能的减少设计变量。在设计优化过程中，减少翼型的几何表示参数可以简化设计过程。因此需要采取合适的翼型参数化方法产生连续光滑的翼型几何外形[3]。

现在典型的参数化方法主要包括Hicks-Henne外形函数法[4]、PARSEC参数化方法[5]、CST翼型参数化方法[6]和样条参数化方法[7]。Hicks-Henne外形函数法具有很强的翼型外形控制能力；PARSEC参数化方法设计参数较少(11个)，具有较好的鲁棒性；CST翼型参数化方法能够描述较大的设计空间；样条参数化方法主要指通过 Bézier曲线、B样条曲线或非均匀有理B样条(NURBS)曲线表示翼型曲线的方法，能够对曲线外形进行局部控制和光滑处理，如 Deng和 Feng[8]在 2011年提出的由4条首尾相接的有理Bézier曲线表示超临界翼型的参数化方法。

随着航空计算的快速发展，上述方法在某些方面已经无法满足现代飞行器外形的精细化设计需求：如Hicks-Henne外形函数法，优化结果存在不光滑现象；PARSEC参数化方法对翼型外形控制能力差，不适用于精细化设计；CST翼型参数化方法鲁棒性差，对超临界翼型拟合效果不理想；而Ferguson样条法[9]所得到的翼型函数不满足C2条件。

B样条方法在几何外形设计中有着广泛的应用[10]，其中三次B样条曲线以其次数不高、满足C2连续性等良好性质而成为最常用的曲线设计方法。利用 B样条方法对于翼型进行参数化设计，曲线的局部调整性有利于之后的翼型优化设计，由于 B样条曲线可以很容易扩展到三维空间，更有利于机翼曲面参数化设计。2008年，Brakhage 和Lamby[11]对B样条方法进行翼型参数化进行了较为详尽的介绍。本文提出一种改进的 B样条翼型参数化方法，在一般的 B样条方法的基础上加入了法向约束，最终用一条四段三次非均匀 B样条曲线表示翼型曲线。得到的翼型曲线的参数个数比一般的 B样条方法有了进一步的减少，所得翼型曲线光滑且满足C2条件。

1 B样条参数化方法

定义 1[12]. 设 n+1个平面向量， Ni,p (t)是定义在节点向量上的p次 B样条基函数(n≥p)，则称：

为相应于节点向量U的p次B样条曲线，称di为控制顶点，顺序连接 d0,…,dn的折线段为控制多边形。

由于三次B样条曲线是应用最广泛的B样条曲线，本文方法采用一条四段三次非均匀 B样条曲线表示翼型曲线。根据如上B样条曲线的定义，曲线 C (t)定义为：

其中定义样条基函数的节点向量为：

其中， u1,u2,u3待定。此节点向量U确保了曲线C (t)在两端插值于首末控制顶点d0与d6。

1.1 目标函数

翼型的参数化就是反求 B样条曲线 C (t)，即根据已知的翼型数据点设计优化目标函数，反求曲线 C (t)的控制顶点，本文采用的优化目标函数为：

表示给定的数据点，C′(ti)表示曲线在点 C (ti)处的切向量，ni表示数据点pi处的单位法向量，X是设计参数向量，即曲线的控制顶点的坐标。根据式(2)，可得：

其中，Nj,3(t)是三次B样条基函数，为C(t)的控制顶点，因此可得：

目标函数 F(X)第一项表示数据点到参数曲线上距离的平方和，保证了参数化曲线与数据点拟合的整体效果；第二项为法向约束，表示数据点的法线向量与参数曲线的切线向量的内积和，确保参数曲线的切线向量与数据点的法向尽量垂直。λ为权系数，可以在区间［0,1］调节，但第一项为主项，λ不宜过大，经过实验将λ取做0.13。

1.2 数据点的法向量估计

为了求解目标函数式(3)，需要首先估计数据点的法向量。简单的估计方法可以用数据点相邻两点的连线的垂线近似数据点的法向，但此方法产生的误差较大。本文利用Jüttler和Felis[13]提出的方法来估计数据点的法向。

如图1所示，方法主要分为2步：

(1) 对于点pi计算局部相关回归线Li，即求解带权的最小二乘问题：

图1 估计数据点法向量

(2) 选取一个新的笛卡尔坐标系，以pi为原点，x轴平行于Li，用二次回归曲线Ci逼近数据点，这条曲线是新坐标系下的一条二次曲线，通过求解：

其中，pnj表示新坐标系下的数据点的坐标，可用二次曲线Ci在原点的法向量近似pi的法向量ni。为了使翼型曲线在前缘满足切线方向垂直于x轴，使最左端前缘数据点的法向平行于x轴，估计出翼型数据点的法向量(图2)。

图2 翼型数据点的法向量

1.3 数据的初步处理

(1) 对给定的翼型数据点进行初步处理，将控制顶点d0设为翼型后缘上表面最后一个数据点的坐标值，再将控制顶点d6设为翼型后缘下表面最后一个数据点的坐标值。剩下需要确定的是控制顶点d1,…,d5。再利用1.2节的方法，估计各数据点对应的法向量。

(2) 由累加弦长法[14]，预估数据点对应的参数，即将数据点对应到参数区间［］0,1上，得到相应参数，并且令t0= 0,tn=1。将节点向量中的u1设为翼型上弧线最高点的参数值，u2设为翼型前缘顶端数据点的参数值，u3设为翼型下弧线最低点的参数值。

1.4 确定控制顶点坐标

通过初步处理已经确定了曲线的首末控制顶点，剩下需要确定的是控制顶点 d1,…,d5，即通过优化1.1节中所给出的目标函数式(3)的过程，使其取得最小值：

其中，等式右端第一项

关于xdk求偏导有(k=1,…,5)：

利用B样条曲线的定义，式(5)右端第二项在节点区间ti∈[Uq,Uq +1]上有(其中Uq表示节点向量第q个节点的值，其中q取3～6)：

其中，

关于xdk求偏导有：

其中

综上，F(X)对关于xdk求偏导有：

同理可得：

令 F(X)关于xdk,ydk(k=1,…,5)的偏导数为零，得到法方程组：

式(6)为含有10个未知数10个方程的方程组，通过求解可以得到控制顶点的坐标从而确定控制顶点。

1.5 误差计算和参数的重新选择

误差是衡量参数化结果好坏的标准，本文中所给的误差均指在翼型弦长为1时各数据点到参数曲线的距离之中的最大值。此外由于数据点的参数化对曲线的形状以及拟合效果具有非常重要的作用，而初始的参数化方法并不一定理想，因此在初始计算的基础上需要对数据点重新参数化并重新计算控制顶点以改善拟合效果，而重新参数化数据点的过程可以与计算拟合误差同步进行。

本文使用基于几何特征的快速迭代法[15]计算点到曲线的距离。如图 3所示，求一点p到曲线C(t)的距离，即希望找到满足 C(t)与p之间距离最小的参数tp，一般根据几何关系可知，矢量ρ=(p-C(tp))必须与曲线在C(tp)点处的切线方向垂直，即满足(p-C (tp))·C′(tp)=0，其中。此方法通过给定初始参数值tp,0，通过迭代逐渐逼近tp。其算法如下：

图3 点到曲线的距离

在计算误差过程中同时更新数据点的参数值，用tp替代数据点的原参数值，再用1.4节的方法求出在新参数值下的新控制顶点，重复1.4节与1.5节的方法，直到误差收敛。此时得到的即为最终的翼型函数。

2 算法与实例

图4为利用三次非均匀B样条曲线进行翼型参数化的算法流程图。采用图 4所示的算法，对EPPLER 398翼型进行参数化。表1给出具体的迭代次数以及每次迭代后的拟合误差。

图4 算法流程图

表1 计算次数与误差

从表1中可以看出，误差在第343次计算达到收敛，至此也就确定了参数曲线的控制顶点坐标。翼型的参数化过程见图5。

图5 EPPLER 398翼型的参数化过程

将这种参数化方法应用于其他一些典型的翼型数据。图6是对NACA 65(2)-215翼型的参数化过程，表 2给出了一些典型翼型的参数化误差，以及参数化的计算时间(算法在 VS2013平台使用C++编程实现，这里的计算时间均是在CPU频率为3.20 Hz，内存为4 GB的PC机上运行得出)。

图6 NACA 65(2)-215翼型的参数化过程

表2 各翼型函数的误差与计算时间

由表1发现，虽然最终翼型函数的误差较小，但误差的收敛速度比较慢。经过实验，对于一般的翼型一般需要迭代 300～400次才能收敛。但本文方法的计算时间是可以控制的，从表 2可以看出，一般翼型的参数化都可以在1 s左右得出结果。

表3 RAE 2822翼型的参数化结果

表4 NACA640-110翼型的参数化结果

3 结论与展望

根据实验结果，对于大多数翼型，本文方法的拟合精度可达到 10–3～10–4，满足机翼外形曲线拟合以及随后生成机翼参数曲面的要求，所得到的翼型函数具有多项优点：光滑且满足C2条件、函数所含的设计变量较少(共 10个)、可以在较短时间内收敛到参数化结果。

但为了达到气动力性能设计的要求，拟合误差需要小于 7×10–4，即翼型的精确表示[11]。将本文方法简单推广到多段(多于四段)[16]，就可以得到一般的改进翼型三次 B样条参数化方法，达到精确表示翼型的要求。表3和表4分别给出了RAE 2822翼型和NACA64-110翼型利用本文算法和典型参数化方法[17]的参数化结果。对比可以发现本文所提出的改进的 B样条参数化方法的参数个数明显少于传统的参数化方法，而将本文方法推广得到的多段改进 B样条参数化方法不仅可以满足气动力性能设计的要求，而且在拟合精度和参数个数方面都比传统的参数化方法有所提高。

在接下来的工作中，将进一步改进算法，增加能量约束[18]，并将算法推广，使得算法可以根据翼型数据点规模和曲率自适应选取 B样条曲线控制顶点个数(即分段数)，从而达到气动力性能设计的要求。最后还将把算法推广到三维上，用 B样条曲面对飞行器翼面进行拟合，从而给出机翼的函数表示。

[1] 刘 周, 朱自强, 付鸿雁, 等. 高升阻比翼型的设计[J]. 空气动力学学报, 2004, 22(4): 410-414.

[2] 王晓璐, 朱自强, 刘 周, 等. 基于N-S方程的翼型双设计点双目标优化设计[J]. 北京航空航天大学学报, 2006, 32(5): 503-507.

[3] Song W, Keane A J. A study of shape parameterization methods for airfoil optim ization, AIAA 2003-4482 [R]. Reston: American Institute of Aeronautics and Astronautics, 2004.

[4] Hicks R M, Hennet P A. Wing design by numerical optim ization [J]. Journal of Aircraft, 1978, 15(7): 407-412.

[5] Sobieczky H. Parametric airfoils and wings [M]//Notes on Numerical Fluid Mechanics (NNFM): Recent Development of Aerodynamic Design Methodologies. Germany: Vieweg Teubner Verlag, 1990: 71-87.

[6] Kulfan B M, Bussoletti J E. “Fundamental” parametric geometry representations for aircraft component shapes, AIAA 2006-6948 [R]. Reston: American Institute of Aeronautics and Astronautics, 2006.

[7] Lepine J, Guibauh F, Trepanier J Y, et al. Optim ized nonuniform rational B-spline geometrical representation for aerodynam ic design of w ings [J]. American Institute of Aeronautics and Astronautics, 2001, 39(11): 570-577.

[8] 邓金秋, 冯仁忠. 利于翼型优化设计的超临界翼型参数化方法[J]. 北京航空航天大学学报, 2011, 37(5): 368-373.

[9] Sobester A, Keane A J. Airfoil design via cubic splines –Ferguson’s curves revisited, AIAA 2007-2881 [R]. Reston: American Institute of Aeronautics and Astronautics, 2007.

[10] 夏成林, 邬弘毅, 郑兴国, 等. 带多个形状参数的三次均匀B样条曲线的扩展[J]. 图学学报, 2011, 32(2): 73-79.

[11] Brakhage K H, Lamby P. Application of B-spline techniques to the modeling of airplane w ings and numerical grid generation [J]. Computer Aided Geometric Design, 2008, 25(9): 738-750.

[12] 王仁宏, 李崇君, 朱春钢. 计算几何教程[M]. 北京:科学出版社, 2008: 158-174, 188-200.

[13] Jüttler B, Felis A. Least-squares fitting of algebraic spline curves via normal vector estimation [J]. Advances in Computational Mathematics, 2002, 17(1): 135-152.

[14] Farin G E. Curves and surfaces for computer-aided geometric design [M]. San Diego: Academ ic Press, 1997: 53-56, 172-224 .

[15] 伍丽峰, 陈岳坪, 谌炎辉, 等. 求点到空间参数曲线最小距离的几种算法[J]. 机械设计与制造, 2011, (9): 15-17.

[16] 韩力文, 杨玉婷, 邱志宇. 一种任意次非均匀 B 样条的细分算法[J]. 图学学报, 2013, 34(5): 56-61.

[17] 张德虎, 席 胜, 田 鼎. 典型翼型参数化方法的翼型外形控制能力评估[J]. 航空工程进展, 2014, 5(3): 281-288.

[18] 徐 进. 带约束的B 样条曲线曲面延伸技术[J]. 图学学报, 2013, 34(3): 36-42.

An Improved Method for Airfoil Parameterization by B-Sp line

Zhang Ji1, Zhu Chungang1, Feng Renzhong2, Liu Mingm ing1, Zhang Hengyang1

(1. School of Mathematical Sciences, Dalian University of Technology, Dalian Liaoning 116024, China; 2. School of Mathematics and Systems Science, Beihang University, Beijing 100083, China)

Airfoil design is a crucial issue of aerodynam ic research, the parameterization of airfoil will affect the airfoil optimization design. In order to reduce the number of variables in the airfoil optim ization, eliminate the unfairness phenomenon, preserve the C2continuity condition, and control the geometric characteristics of the airfoil in the optim ization process, in this paper, we present an improved method for airfoil parameterization by B-spline. The method represents airfoil by a cubic non-uniform B-spline curve. Fitting of airfoil data by B-spline curve is mainly by least square method and the normal constraints. And the final result is obtained by iteration. Experiments show that the proposed method can be well fitted to the typical airfoil data, the resulting curve is fair and C2continuity, and has few parameters of airfoil function compared with the classical airfoil parametric methods.

airfoil; B-spline curves; airfoil parameterization; airfoil optim ization

TP 391

10.11996/JG.j.2095-302X.2016030342

A

2095-302X(2016)03-0342-07

2015-11-05；定稿日期：2015-12-01

国家自然科学基金项目(11271060, 11290143)；民用飞机专项项目(MJ-F-2012-04)；中央基本科研业务费资助项目(DUT16LK38)；辽宁省高等学校优秀人才支持计划项目(LJQ2014010)

张 骥(1990–)，男，江苏盐城人，硕士研究生。主要研究方向为计算几何。E-mail：jizhang1990@gmail.com

朱春钢(1977–)，男，北京人，教授，博士。主要研究方向为计算几何与计算机辅助几何设计。E-mail：cgzhu@dlut.edu.cn