2. 中国空气动力研究与发展中心, 四川 绵阳 621000
2. China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China
0 引 言
随着计算机以及CFD[1, 2]技术的高速发展，用数值方法进行飞行器的气动外形优化设计[3]已成一种趋势。在气动外形优化设计中，外形的参数化是非常重要的一个环节，好的参数化方法不仅能提高优化效率而且能得到很好的结果。在翼型参数化的过程中，主要关心设计变量的个数以及设计空间的大小，少的设计变量有助于提高优化算法的收敛速度，充足的设计空间有助于找到可能的最优解，但二者互相矛盾。直接使用较少的设计变量可以提高优化搜索效率，但是可能由于设计空间小、设计几何外形变化生硬、缺乏必要的灵活性等，导致难以获得优良的结果；直接使用多的设计变量虽然设计空间充足，但是由于各个设计变量随机盲目变化，几何外形变化的过于灵活，中间过程可能会产生很多畸形无用的外形，导致优化搜索的效率大大降低，甚至难以获得收敛的优化设计结果。
本文以翼型设计为例，讨论了不同阶次Bezier样条曲线对超临界翼型的描述能力，然后基于该曲线发展了一套逐次升阶的翼型参数化方法，结合改进的粒子群优化算法，构建了逐步扩展设计空间的气动优化设计方法。该方法首先利用低阶Bezier曲线参数化构建的设计空间，优化搜索出性能较为优秀的翼型。在此基础上，利用曲线的升阶保形特性，对曲线进行升阶，传递设计变量，扩展更新粒子群，继续进行优化搜索，不断地增加设计变量个数，不断地扩大设计空间。由于设计变量个数由少到多，低阶优化搜索结果为高阶曲线指明了“方向”，使得该优化设计方法在提高优化质量的同时，提高了优化搜索的效率。 1 基于低阶Bezier曲线的翼型参数化
Bezier曲线[4, 5]是一种参数样条曲线。不同于插值样条曲线，它是一种逼近多边形顶点的方法。Bezier曲线用参数方法表示自由曲线，具有几何不变性，可以处理无穷大斜率和多值曲线，并易于进行坐标变换等诸多优点。用Bernstein多项式作为基函数描述Bezier曲线的表达式为：
式中Vi(i=0,1,…,n)为特征多边形的顶点，Bn,i(u)=Cni(1－u)n－iui(i=0,1,…,n)为Bernstein多项式，u∈[0,1]为变化参数。顶点的个数为n+1，而且顶点的位置一旦确定，曲线的形状就唯一确定。
Bernstein多项式具有非负性、混合性，
以及对称性，Bn,i(u)=Bn,n－i(1－u)等特点。因此，它在自由曲线、曲面建模中得到广泛的应用。
由于超临界翼型下表面通常存在拐点，所以能很好地考验翼型参数化方法的能力。本文首先采用低阶(n=5)Bernstein多项式，对RAE2822超临界翼型采用最小二乘法进行拟合。翼型上下表面各采用5次Bezier曲线，各 6个控制点(其中头部和尾部的顶点保持不动)，设计变量为顶点的纵坐标，共8个(上下表面各4个设计变量)。
图 1显示了五次Bernstein多项式的函数图。由Bernstein多项式的混合性，可知：对任意给定参数u，Bezier曲线就是特征多边形顶点Vi的加权平均，其权值因子就是Bernstein多项式的值。
 |
| 图 1 五次Bernstein多项式函数图Fig. 1 Figure of fifth-order Bernstein polynomial |
图 2为五次Bezier曲线拟合RAE2822翼型充分收敛[6]后的几何外形。图 3给出了上下翼型表面的拟合误差。该算例中使用8个设计变量的Bezier曲线，描述翼型的最大几何误差为0.0011。为了检验翼型外形偏差对其气动特性的影响，本文利用CFD方法对拟合翼型的流场与气动特性进行了计算，并将计算结果与原始翼型RAE2822的计算结果进行了比较。计算状态为：
Ma=0.5,α=2°,Re=5×106
 |
| 图 2 五次Bezier曲线拟合RAE2822Fig. 2 Approximation of RAE2822 using fifth-order Bezier curves |
 |
| 图 3 上下表面拟合误差Fig. 3 Approximation residuals of the upper and lower surface |
表 1为拟合翼型与RAE2822翼型的几何误差与气动性能偏差。由表 1可见，两个翼型的气动性能存在偏差，升力系数、阻力系数和力矩系数相对误差的平均值为2.24%。假设在某次优化设计中，RAE2822就是需要的“最优”翼型，8个设计变量的Bezier曲线由于变化空间有限，就不能很好的描述这个“最优”翼型。换言之，为了获得最佳的外形，必须提高曲线的参数化能力，扩大设计空间进而提高对某个翼型描述的精度。
| RAE2822 | Fitted airfoil | Relative error | |
| CL | 0.483 | 0.497 | 2.90% |
| CD | 0.0095 | 0.0094 | 0.37% |
| Cm | 0.072 | 0.074 | 3.45% |
| Average | 2.24% |
通过对RAE2822的参数化拟合，发现用较少的设计变量对翼型参数化时，每个设计变量的变化对曲线形状的影响程度过大，对曲线描述不够灵活，不能精准地描述翼型。因此，当用于翼型设计时，由于设计空间有限，导致设计结果质量不高。
文献[6]提出了广义Bezier曲线[7, 8]的概念，即在给定控制顶点和相应Bernstein基函数的情况下，曲线的形状仍然可以再微调。图 4显示了几组广义Bezier曲线，该组曲线只是单纯扩大了曲线的变化范围，但它在精确描述某个翼型方面能力有限。
 |
| 图 4 广义Bezier曲线Fig. 4 General Bezier curves |
为解决此问题，文中尝试曲线升阶的方法，即在保持原低阶拟合曲线形状不变的情况下，增加控制点的数量，提高曲线的阶次，增加曲线的柔性和变化的灵活性。
基于Bernstein基函数的升阶算法表示为：
上式中n为曲线在升阶前的幂次，Vni为升阶前的控制点。此法应用于递推公式，每次增加一个控制点(升一阶)，图 5显示了二次曲线升到三次曲线的情形。由文献[4]可知，低阶的Bezier曲线一定能用高阶的Bezier曲线来表示，反之则不可以，即高阶曲线的设计空间一定大于低阶曲线。
 |
| 图 5 二次曲线升阶到三次曲线示意图Fig. 5 Sketch map of two-order curve and three-order curve |
在低阶曲线(8个设计变量)拟合结果的基础上，对曲线进行升阶，分别逐步升到12、16、20和24个设计变量，如图 6~图 9所示。
 |
| 图 6 8个设计变量升到12个设计变量Fig. 6 Increasing variables from 8 to 12 |
 |
| 图 7 8个设计变量升到16个设计变量Fig. 7 Increasing variables from 8 to 16 |
 |
| 图 8 8个设计变量升到20个设计变量Fig. 8 Increasing variables from 8 to 20 |
 |
| 图 9 8个设计变量升到24个设计变量Fig. 9 Increasing variables from 8 to 24 |
图 10显示了升阶后不同设计变量拟合RAE2822翼型的误差，发现曲线阶次的增加明显提高了曲线拟合的几何精度。图 11和图 12分别显示了不同设计变量收敛后翼型上下表面的拟合误差。 图 13和图 14分别显示了不同设计变量收敛后翼型上下表面曲线的曲率半径归1化后的分布。由此可见，随着设计变量的4增加，拟合误差越来越小，但是翼型表面曲率的波动越来越明显。
 |
| 图 10 不同设计变量的拟合误差Fig. 10 Approximation error of different variables |
 |
| 图 11 不同设计变量拟合后上表面误差Fig. 11 Approximation error of the upper surface of different variables |
 |
| 图 12 不同设计变量拟合后下表面误差Fig. 12 Approximation error of the lower surface of different variables |
 |
| 图 13 翼型上表面曲率半径分布Fig. 13 Radius of curvature distribution of upper surface |
 |
| 图 14 翼型下表面曲率半径分布Fig. 14 Radius of curvature distribution of lower surface |
表 2对比了不同设计变量下的气动性能和相对误差的平均值(与上文的飞行状态相同)，发现超过20个(包括20)设计变量时，气动力系数的平均相对误差小于1%，达到了很高的精度。同时发现，设计变量达到24个时，由于曲线拐点增多，波动性加强，很难满足翼型表面光顺性的要求，气动误差较20个设计变量时反而变大。对于翼型优化来说，综合考虑拟合误差和曲线的曲率分布，最终使用20个设计变量足够。
| CL | CD | Cm | Average of relative error | |
| RAE2822 | 0.483 | 0.0095 | 0.072 | |
| 8 variables | 0.497 | 0.0094 | 0.074 | 2.24% |
| 12 variables | 0.470 | 0.0093 | 0.068 | 3.11% |
| 16 variables | 0.473 | 0.0094 | 0.069 | 1.99% |
| 20 variables | 0.480 | 0.0094 | 0.071 | 0.65% |
| 24 variables | 0.487 | 0.0094 | 0.072 | 0.81% |
基于Bezier曲线的特性，建立逐步升阶的翼型参数化方法：在优化搜索的过程中，给曲线逐次升阶，设计变量由少到多，不断地扩大设计空间。 3 改进的粒子群算法及优化设计方法
在粒子群算法[9, 10](PSO)的可调参数中，惯性权值是最重要的参数，较大的w有利于提高算法的全局搜索能力，而较小的w会增强算法的局部搜索能力。文中采用改进后的自适应惯性权值[11, 12]粒子群算法(SAPSO)进行优化搜索。自适应权重法的惯性权重系数公式表达式如下：
其中wmax、wmin分别表示w的最大值和最小值，f表示当前的目标函数值，fave、fmin分别表示当前所有粒子的平均目标值和最小目标值。式(3)中，目标函数值优于平均目标值的微粒，其对应的惯性权重因子较小，从而保护该微粒；反之，对于目标函数值差于平均目标值的微粒，其对应的惯性权重因子较大，使得该微粒去寻找较好的搜索区域。
基于改进的粒子群算法，若直接使用20个设计变量进行优化搜索，设计变量过多且变化随机盲目，导致搜索效率大大降低。本文结合逐步升阶的翼型参数化方法和改进的粒子群算法，构建逐步扩展设计空间的优化设计方法。首先在低阶曲线参数化构建的设计空间中，进行粒子群优化搜索，进化一定代数后，保持已优化外形不变的情况下，按照上文的算法给曲线升阶，增加设计变量的个数，进行设计变量的传递，扩展更新粒子群(即增加每个粒子的维度)，继续进行粒子群搜索。即先由较低阶曲线寻找到有利的区域，然后增加设计变量个数，用高阶曲线在低阶的基础上继续在该区域精细化寻优，逐层优化搜索。这样，每次升阶的过程都保留低阶的最优结果，不但扩大设计空间，还能提高优化搜索的效率。图 15给出了该优化设计方法的框架。
 |
| 图 15 优化设计方法框架Fig. 15 Frame of the optimization method |
本文的气动分析方法采用基于结构网格的RANS方程[13, 14]，网格数量为4万，使用Roe空间离散方法，S-A一方程湍流模型，LU-SGS隐式时间推进方法，多重网格加速收敛。为提高计算机的并行效率，对结构网格进行分块处理。
以NASA0412超临界翼型为初始翼型进行减阻优化设计，设计状态为：
Ma=0.75,CL=0.69,Re=1×107
优化目标为最小化翼型的阻力系数，约束条件为厚度不减小，力矩系数绝对值不大于0.12，不满足约束则对目标函数进行惩罚。优化的数学模型可描述为：
为了验证本文方法的优越性，对比本文优化方法与传统Hicks-Henne型函数翼型参数化的优化方法如下：
1) 参数设置。本文方法：初始为8个设计变量，每次升阶增加4个设计变量，最后为20个设计变量；Hicks-Henne方法：使用20个设计变量，设计变量个数不变。
2) 优化算法设置。本文方法：每代粒子群数目为50，不同的设计变量各进化10代，共40代；Hicks-Henne方法：每代粒子群数目50，连续进化40代。
图 16和图 17分别给出了本文方法优化前后的翼型形状和翼型表面的压力分布。在设计状态下，初始翼型上面出现较强的激波，引起很大的激波阻力。优化后的翼型下表面前缘半径减小，压力增加，前加载增大；最大厚度点后移，弯度增加，后加载增大；由于升力系数固定，前后加载增大，翼型上表面的压力缓慢恢复，激波基本消除，大大减小了翼型的阻力。
 |
| 图 16 优化前后的翼型形状Fig. 16 Comparison between initial and optimized airfoil configuration |
 |
| 图 17 优化前后的压力分布Fig. 17 Comparison of pressure distribution between initial and optimized airfoil |
表 3给出了使用本文方法和使用Hicks-Henne方法的优化前后的结果对比，本文的方法优化阻力系数减小了0.0057，优化质量高于Hicks-Henne方法的0.0054。
| Initial | Opt.(adding) | Opt.(Hicks-Henne) | |
| CL | 0.69 | 0.69 | 0.69 |
| CD | 0.0179 | 0.0122 | 0.0125 |
| Cm | -0.105 | -0.106 | -0.108 |
| Thick | 0.12 | 0.1206 | 0.1205 |
图 18给出了两种方法优化收敛历程对比，很明显文中的翼型优化方法在搜索效率和质量方面高于传统的Hicks-Henne参数化的优化方法。使用文中的参数化方法，每次设计变量的增加时(每10代增加一次)，优化的质量都提高一次，打破原来的收敛趋势，证明了该参数化方法的有效性和实用性。
 |
| 图 18 优化收敛历程对比Fig. 18 Comparison of convergence history |
文中研究了翼型参数化中设计变量和设计空间关系，基于Bernstein多项式的性质，提出Bezier曲线逐次升阶的翼型参数化方法，结合改进的粒子群算法构建了逐步扩展设计空间的优化设计系统，很好地解决了设计变量和设计空间这对矛盾。最后，通过对典型超临界翼型的减阻设计，对比了文中方法和传统的优化方法，证实了该方法的可行性和高效性。同时，该升阶思想可以拓展到三维Bezier曲面建模以及以Bernstein多项式为基函数的自由变形造型(FFD)[15, 16]中，有很好的实用价值。
| [1] | John D Anderson. Computational fluid dynamics[M]. Beijing: Tsinghua University Press, 2005. |
| [2] | Zhu Ziqiang. Application of computational fluid dynamics[M]. Beijing: Beihang University Press, 1998. (in Chinese) 朱自强. 应用计算流体力学[M]. 北京: 北京航空航天大学出版社, 1998. |
| [3] | Su Wei. Aerodynamic optimization design based on computationalfluid dynamics and surrogate model[D]. Xi'an: Northwestern Polytechnical University, 2007. (in Chinese) 苏伟. 基于CFD技术和代理模型的气动外形优化设计方法研究[D]. 西安: 西北工业大学, 2007. |
| [4] | Zhu Xinxiong. Free curve and surface modeling techniques[M]. Beijing: Science Press, 2000. (in Chinese) 朱心雄. 自由曲线曲面造型技术[M]. 北京: 科学出版社, 2000. |
| [5] | Zhu Xinxiong. Lecture notes on representations of curves and surfaces[M]. Minnesota: University of Minnesota, 1981. |
| [6] | Jin Yaochu. A comprehensive survey of fitness approximation in evolutionary computation[J]. Soft Computing, 2005, (9): 3-12. |
| [7] | Lu Yueqi. General Bezier curves with alterable free parameters[D]. Hangzhou: Zhejiang University, 2006. (in Chinese) 卢跃奇. 具有可调自由参数的广义Bezier曲线[D]. 杭州: 浙江大学, 2006. |
| [8] | Zhu Xiumei. Research on some parameter curves[D]. Hefei: HefeiUniversity of Technology, 2006. 朱秀梅. 几类参数曲线的研究[D]. 合肥: 合肥工业大学, 2006. |
| [9] | Wang Rongwei. Research on optimization algorithm and application in aerodynamic optimization design[D]. Xi'an: Northwestern Polytechnical University, 2011. (in Chinese) 王荣伟. 优化算法及其在气动优化中的应用研究[D]. 西安: 西北工业大学, 2011. |
| [10] | Gong Chun, Wang Zhengling. Master MATLAB for optimization [M]. Beijing: Publishing House of Electronics Industry, 2009. (in Chinese) 龚纯, 王正林. 精通MATLAB最优化设计[M]. 北京: 电子工业出版社, 2009. |
| [11] | Li Mingwei. Partical swarm optimization algorithm and application in aerodynamic optimization design[D]. Xi'an: Northwestern Polytechnical University, 2009. (in Chinese) 黎明伟. 粒子群优化算法研究及其在气动力优化设计中的应用[D]. 西安: 西北工业大学, 2009. |
| [12] | Li Ding. Application of improved partical swarm optimization algorithm to aerodynamic design[J]. Acta Aeronoutica et Astronautica Sinica, 2012, 33(10): 1809-1816. (in Chinese) 李丁. 改进的粒子群优化算法在气动设计中的应用[J]. 航空学报, 2012, 33(10): 1809-1816. |
| [13] | Huang Jiangtao. Research on aircraft aerodynamic configuration optimization design method and applications[D]. Xi'an: Northwestern Polytechnical University, 2012. (in Chinese) 黄江涛. 飞行器气动外形优化设计方法及应用研究[D]. 西安: 西北工业大学, 2012. |
| [14] | Li Jing. Aerodynamic optimization system based on CST technique[J]. Acta Aerodynamica Sinica, 2012, 30(4): 443-449. (in Chinese) 李静. 基于CST参数化方法气动优化设计研究[J]. 空气动力学学报, 2012, 30(4): 443-449. |
| [15] | Huang Jiangtao. Laminar airfoil aerodynamic optimization design based on delaunay graph mapping and FFD technique[J]. Acta Aeronoutica et Astronautica Sinica, 2012, 33(10): 1817-1826. (in Chinese) 黄江涛. 应用Delaunay图映射与FFD技术的层流翼型气动优化设计[J]. 航空学报, 2012, 33(10): 1817-1826. |
| [16] | Jamshid A. Aerodynamic shape optimization based on free-form deformation[R]. AIAA 2004-4630. |