基于FFD方法的高超声速升力体气动优化

夏陈超，邵 纯，姜婷婷，陈伟芳

(浙江大学 航空航天学院,杭州 310027)



基于FFD方法的高超声速升力体气动优化

夏陈超，邵 纯，姜婷婷，陈伟芳

(浙江大学 航空航天学院,杭州 310027)

将FFD(Free Form Deformation)自由变形法与无限插值动网格方法相结合，发展了一种飞行器参数化建模和网格生成方法。二维和三维的实例显示自由变形之后得到的飞行器几何外形及其对应的网格能保持平滑光顺，验证了方法的有效性。在此基础上，结合径向基函数代理模型和CFD技术发展了一套优化设计方法并对高超声速升力体外形进行了气动优化。基于自适应模拟退火算法的单目标优化表明，在保持原有外形体积不减小的情况下，升阻比提高了1.28%；基于NSGA-II的多目标优化得到了飞行器升阻比和体积的最优解集，典型优化外形的升阻比和体积分别提高了2.93%和2.49%。升力体的优化结果表明了FFD方法的有效性和优化设计方法的实用性。

自由变形；无限插值法；高超声速；升力体；优化设计

0 引言

随着计算流体力学(CFD, Computational Fluid Dynamics)技术的发展和计算机水平的提高，基于CFD的高保真度优化在飞行器设计中日渐增多。在优化设计过程中，参数化建模和网格生成是两个基础而又十分重要的环节。高效的优化设计系统应该包含快速、便捷和适用性强的参数化建模和网格生成模块。常用的飞行器参数化建模方法有基于CAD法、解析法、多项式/样条曲线法、自由变形法(FFD, Free Form Deformation)等[1]。相比于传统的参数化方法，FFD方法可用较少的参数实现对复杂几何外形的控制，操作便捷，变形能力强，同时能够保持原有外形的光滑和连续，是较为理想的参数化方法。Taira等[2]基于FFD技术，建立了可视化的飞行器布局设计平台；Lyu等[3]采用FFD方法对翼身融合体进行了极多设计变量的气动外形优化，降低了飞行器阻力。国内，王丹等[4]、黄江涛等[5]和王元元等[6]采用FFD技术分别对翼身组合体整流包、翼稍小翼和运输机后体进行了优化设计，取得了良好的效果。除了参数化建模，为提高高保真度优化设计系统的效率，需要在几何外形发生变化时，迅速形成新的网格用于数值计算。新的网格可以在优化过程中重新生成[7]，更多的是采用动网格技术生成，原因是复杂外形飞行器的网格重新自动生成十分困难和耗时。常用的动网格方法有无限插值法[8-9]、弹簧法[10]和Delaunay图映射法[11]等，其中无限插值法是一种简单高效的结构动网格技术，在飞行器优化设计、非定常数值计算中得到了广泛应用。升力体是典型的高超声速滑翔飞行器外形，具有较好的升阻特性，Kinney[12]和张珍铭等[13]采用工程估算方法分别对HL-20和类X-33升力体进行了优化设计，马洋等[14]采用CFD方法进行了滑翔式升力体外形优化，均在一定程度上提高了飞行器升阻比。在采用CFD方法进行飞行器优化时，计算量往往十分巨大，一种折中措施是采用代理模型的近似方法，Kriging模型[15]、径向基函数模型(RBF, Radial Basis Function)[16]和响应面模型(RSM, Response Surface Model)[17]等是常用的代理模型，在飞行器设计中应用较多。

目前FFD方法主要应用于亚声速流动，在高超声速流动中的应用还很少。本文将FFD自由变形法与无限插值动网格方法相结合，建立了一种参数化建模和网格生成方法，接着基于径向基函数代理模型、CFD方法和单目标/多目标优化算法，对高超声速升力体进行了气动外形优化设计。

1 FFD自由变形和动网格方法

FFD自由变形法由Sederberg和Parry[18]于1986年提出，其基本思想是将几何模型嵌入到一个控制体中，建立数学关系后通过改变控制顶点来改变模型的形状。本文采用Bernstein多项式来定义控制体，若在原点X0建立一个局部坐标系，则控制体中任意一点的全局坐标可写为

X=X0+sS+tT+uU

(1)

其中，s、t、u为该点在局部坐标系中沿S、T、U方向上的值，可表示为

(2)

若局部坐标系3个轴矢量方向上分别有l、m、n个控制顶点，则控制体将含(l+1)(m+1)(n+1)个顶点，其全局坐标可表示为

(i=0,1…,l；j=0,1,…,m；k=0,1,…,n)

(3)

对于待变形模型上任意一点，在求得其局部坐标(s,t,u)后，变形后新的全局坐标可由Bernstein多项式计算得到

(4)

实施FFD自由变形的主要步骤为：根据待变形几何模型建立控制体和局部坐标系；计算几何模型的局部坐标；改变控制顶点的全局坐标；根据式(4)计算几何模型新的全局坐标。

在完成几何变形之后，即可采用TFI方法进行计算网格的变形。TFI无限插值法[19]基于原有的多块对接结构网格，通过归一化的弧长参数将变形量传递至流场网格。以三维的体网格插值为例，其无限插值过程为

ΔSi,j,k=U+V+W-UV-UW-VW+UVW

(5)

其中

U=(1-αi,j,k)ΔS(1,j,k)+αi,j,kΔS(NI,j,k)

V=(1-βi,j,k)ΔS(i,1,k)+βi,j,kΔS(i,NJ,k)

W=(1-γi,j,k)ΔS(i,j,1)+γi,j,kΔS(i,j,NK)

UV= (1-αi,j,k)(1-βi,j,k)ΔS(1,1,k)+

αi,j,k(1-βi,j,k)ΔS(NI,1,k)+

(1-αi,j,k)βi,j,kΔS(1,NJ,k)+

αi,j,kβi,j,kΔS(NI,NJ,k)

UW= (1-αi,j,k)(1-γi,j,k)ΔS(1,j,1)+

αi,j,k(1-γi,j,k)ΔS(NI,j,1)+

(1-αi,j,k)βi,j,kΔS(1,j,NK)+

αi,j,kβi,j,kΔS(NI,j,NK)

VW= (1-βi,j,k)(1-γi,j,k)ΔS(i,1,1)+

βi,j,k(1-γi,j,k)ΔS(i,NJ,1)+

(1-βi,j,k)γi,j,kΔS(i,1,NK)+

βi,j,kγi,j,kΔS(i,NJ,NK)

UVW= (1-αi,j,k)(1-βi,j,k)(1-γi,j,k)ΔS(1,1,1)+

αi,j,k(1-βi,j,k)(1-γi,j,k)ΔS(NI,1,1)+

(1-αi,j,k)βi,j,k(1-γi,j,k)ΔS(1,NJ,1)+

(1-αi,j,k)(1-βi,j,k)γi,j,kΔS(1,1,NK)+

(1-αi,j,k)βi,j,kγi,j,kΔS(1,NJ,NK)+

αi,j,k(1-βi,j,k)γi,j,kΔS(NI,1,NK)+

αi,j,kβi,j,k(1-γi,j,k)ΔS(NI,NJ,1)+

αi,j,kβi,j,kγi,j,kΔS(NI,NJ,NK)

式中 ΔS为位移；α、β、γ为归一化的弧长参数；NI、NJ、NK为网格块3个方向编号的最大值。

2 自由变形和动网格实例

2.1 二维RAE2822翼型

RAE2822是经典的二维跨声速翼型，常作为湍流模型的验证算例。图1为原始的翼型网格和FFD控制体，网格在近壁面进行了加密，圆点为控制顶点。图2为控制顶点、翼型的变化以及对应的网格。可看出，翼型和网格在变形之后均可保持连续和光顺。

图1 原始网格和控制体Fig.1 Initial grid and control volume

图2 变形后网格和控制体Fig.2 Grid and control volume after deformation

2.2 三维NACA64A010机翼

以三维NACA64A010机翼[20]进行验证，图3为机翼翼尖扭转15°后的网格和拓扑，图4为机翼向上翘曲的网格示意图，其中圆点为控制顶点，控制体内为变化后外形。可以看出，FFD方法的变形能力较强，且变形后的机翼网格仍能保证较高的质量。

图3 翼尖扭转15°前后的网格Fig.3 Comparison of grid before and after 15°twist of wing tip

图4 机翼翘曲前后的网格Fig.4 Comparison of grid before and after the warp of wing

2.3 类X-33升力体外形

以类X-33升力体[21]外形为例，检验自由变形和动网格方法对复杂飞行器外形的适用性。图5为飞行器半模网格，网格量约140万，共含29个网格块。图6为飞行器头部迎风面局部变形示意图，控制顶点的移动改变了飞行器表面局部形状。可以看出，飞行器几何外形在经过FFD变形后仍能保持很好的光滑性。FFD技术的这一特点使其十分适合于飞行器的局部修形和精细优化。

图5 类X-33网格Fig.5 Grid of the X-33-type vehicle

图6 类X-33头部迎风面FFD变形Fig.6 Free form deformation of the windward side of the X-33-type vehicle

3 升力体外形优化设计

基于发展的自由变形和动网格技术，构建优化系统并针对一类简化的高超声速升力体外形进行优化设计。升力体初始外形长2 m，底部为长、短半轴之比为3∶1的椭圆形。初始网格如图7所示，网格量约60万，划分为18块，网格在壁面进行了加密保证y+值小于1。包围升力体的FFD控制体如图8所示，放大图为控制底部截面形状的控制顶点及其编号。设计变量为1～6号控制点(考虑了左右对称性)，保持飞行器长宽不变，控制点仅在Y方向移动。

图7 升力体网格示意图Fig.7 Grid of the lifting body

图8 控制体和设计变量示意图Fig.8 Control volume and design variables

3.1 数值计算方法

飞行器气动特性的获取通过求解可压缩N-S方程，对流项采用二阶TVD格式进行离散，时间项采用隐式推进方法，湍流模型为k-εSST[22]两方程模型。进口为来流条件，出口为外推条件，壁面为等温条件。

3.2 试验设计和代理模型

构建合适的代理模型是解决高精度优化中计算量过大的常用措施，其过程包括试验设计(DOE, Design of Experiment)、样本点计算以及代理模型的建立和验证。本文升力体的设计变量如表1所示。基于此，采用拉丁超立方(LHS, Latin Hypercube Sampling)试验设计方法形成40个样本点，结合适用于非线性全局优化问题的径向基函数方法构造代理模型。

表1 设计变量范围Table 1 Design variables and ranges

图9给出了飞行器体积和升阻比的分析值与预测值之间的对比。其中，体积的分析值通过飞行器体积分得到，升阻比的分析值通过CFD计算得到，预测值均通过代理模型获得。通过交叉验证方式计算得到代理模型对飞行器体积和升阻比的均方根误差(RMSE, Root Mean Square Error)分别为2.6%和3.3%，决定系数(R-Squared)分别为98.5%和98%。可以看出，构建的代理模型具有较高的精度，满足优化设计的需求。

(a)体积 (b)升阻比

3.3 优化流程和算法

本文基于代理模型方法，构建了优化设计系统，基本流程如图10所示。在建立满足精度要求的代理模型后，采用自适应模拟退火算法(ASA, Adaptive Stimulated Annealing)[23]进行升阻比的单目标优化，采用非支配排序遗传算法(NSGA-II, Non-dominated Sorting Genetic Algorithm)[24]进行升阻比和体积的多目标优化。

图10 优化设计流程Fig.10 Flowchart of optimization

3.4 优化结果与分析

针对本文的升力体外形，在马赫数10，攻角10°，高度30 km条件下进行优化，初始外形升阻比2.186，体积0.281m3。首先进行单目标优化，约束条件为飞行器体积不小于初始外形。图11为采用自适应模拟退火算法优化后升阻比的收敛过程，最优外形对应的升阻比为2.214，相比于初始外形提高了1.28%。通过对最优外形进行CFD计算得到的升阻比为2.203，与预测值相对误差仅为0.5%。从图12的底部截面形状对比可看出，最优外形在飞行器对称面上的厚度略有增加，在对称面与翼前缘之间的厚度有一定减小。

用NSGA-II优化算法进行体积和升阻比的多目标优化，种群数量为400，进化代数为100。由于体积和升阻比相互制约，因而多目标优化得到的是关于体积和升阻比的诸多可行解及最优解集(Pareto前沿)，如图13。图中三角形点为初始外形，空心原点为Pareto前沿上的4个典型结果。对于外形1结果，其体积为0.676 m3，但升阻比只有0.95，对于外形4结果则相反，其升阻比虽高达2.83，但体积只有0.103 m3。

图13中外形3的升阻比和体积分别为2.250和0.288 m3，采用CFD数值计算得到的升阻比和体积分别为2.199、0.290 m3，误差分别仅为2.32%和0.69%。相比于初始外形的升阻比和体积分别提高了2.93%和2.49%。图14为外形3底部截面与初始外形的对比，与单目标优化结果类似，此时飞行器迎风面的对称面位置明显增厚，翼前缘与对称面之间减薄，即图8中序号为2的设计变量向下移动量较大，序号为6的设计变量向上移动量较大。从图15可看出外形3升阻比比初始外形大的原因：外形3翼前缘附近的马赫数等值线(实线)比初始外形(虚线)更加贴近前缘，即其下表面的高压气流更少地溢流到上表面，有类似于乘波体的效果，因而升阻比较大。综合单目标和多目标优化，对本文10°攻角的优化状态，迎风面形状对升阻比影响很大。实际飞行器设计中应根据需求，按照不同权重对升阻比和体积折中选择。

图11 升阻比收敛过程Fig.11 Convergence history of lift-drag ratio

图12 优化外形和初始外形的对比Fig.12 Comparison of initial and optimized shapes

图13 多目标Pareto前沿Fig.13 Pareto front of multi-objectiveo ptimization

图14 外形3和初始外形的对比Fig.14 Comparison of initial and optimized shapes of configuration 3

图15 翼前缘附近马赫数等值线对比Fig.15 Comparison of Mach contours near leading edge

4 结论

(1)基于FFD自由变形法和TFI无限插值技术，发展了参数化建模和网格生成方法，通过对翼型、机翼和复杂外形飞行器的建模与动网格生成验证了方法的有效性。FFD方法操作便捷，变形能力强，能保持几何外形的光顺，TFI方法简易高效，二者的结合可为优化设计提供良好的基础，提高优化效率。

(2)基于径向基函数代理模型和CFD技术，发展了高保真度优化设计方法，针对高超声速升力体外形进行了单目标/多目标优化设计。在保持初始外形体积不变的条件下，基于ASA算法的单目标优化将升阻比提高了1.28%。采用NSGA-II算法得到了升力体升阻比和体积的最优解前沿，最优解典型外形的升阻比和体积较初始外形分别提高了2.93%和2.49%，相比于初始外形其迎风面腹部厚度增大，腹部两侧厚度减小。

(3)本文发展的优化设计方法基本满足飞行器优化设计的需求，具有较好的工程实用性。在此基础上将进一步开展飞行器气动布局优化设计工作。

[1] Samareh J A.Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization[J].AIAA Journal,2001,39(5):877-884.

[2] Taira T,Jeong S,Obayashi S,et al.GUI-based geometry deformation tool for modification of aircraft configurations[R].AIAA 2013-0619.

[3] Lyu Z,Martins J.RANS-based aerodynamic shape optimization of a blended-wing-body aircraft[R].AIAA 2013-2586.

[4] 王丹,白俊强,黄江涛.FFD方法在气动优化设计中的应用[J].中国科学:物理学力学天文学,2014,44(3):267-277.

[5] 黄江涛,高正红,白俊强,等.基于任意空间属性FFD技术的融合式翼稍小翼稳健型气动优化设计[J].航空学报,2013,34(1):37-45.

[6] 王元元,张彬乾,郭兆电,等.基于FFD技术的大型运输机上翘后体气动优化设计[J].航空学报,2013,34(8):1806-1814.

[7] 夏陈超,赵文文,陈伟芳,等.类HTV-2升力体参数化建模与网格自动生成研究[C]//第十五届全国激波与激波管学术交流会文集.杭州,2012:503-508.

[8] Byun C,Guruswamy G P.A parallel,multi-block,moving grid method for aeroelastic applications on full aircraft[R].AIAA 98-4782.

[9] Tsai H M,F.Wong A S,Cai J,et al.Unsteady flow calculations with a parallel multiblock moving mesh algorithm[J].AIAA Journal，2001,39(6):1021-1029.

[10] Zeng D,Ethier C R.A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains[J].Finite Elements in Analysis and Design,2005,41(11):1118-1139.

[11] Liu X,Qin N,Xia H.Fast dynamic grid deformation based on Delaunay graph mapping[J].Journal of Computational Physics,2006,211(2):405-423.

[12] Kinney D.Aerodynamic shape optimization of hypersonic vehicles[R].AIAA 2006-239.

[13] 张珍铭,丁运亮,刘毅.升力体外形设计的代理模型优化方法[J].宇航学报,2011,32(7):1435-1444.

[14] 马洋,杨涛,张青斌.高超声速滑翔式升力体外形设计与优化[J].国防科技大学学报,2014,36(2):34-40.

[15] 杨希祥,周张,彭科.基于Kriging方法的整流罩气动外形设计优化[J].固体火箭技术， 2014,37(2):167-171.

[16] Leary S J,Bhaskar A,Keane A J.Global approximation and optimization using adjoint computational fluid dynamics codes[J].AIAA Journal,2004,42(3):631-641.

[17] 成沉,鲍福廷,刘旸,等.基于响应面法的喉栓式喷管型面优化设计[J].固体火箭技术,2014,37(1):47-51.

[18] Sederberg T W,Parry S R.Free-form deformation of solid geometric models[J].Computer Graphics,1986,20(4):151-160.

[19] Thompson J F,Soni B K,Weatherill N P.Handbook of grid generation[M].CRC Press,1998.

[20] Chaderjian N M,Guruswamy G P.Transonic Navier-Stokes computations for an oscillating wing using zonal grids[J].Journal of Aircraft,1992,29(3):326-335.

[21] Murphy K J,Nowak R J,Thompson R A,et al.X-33 hypersonic aerodynamic characteristics[J].Journal of Spacecraft and Rockets,2001,38(5):670-683.

[22] Menter F R,Kuntz M,Langtry R.Ten years of industrial experience with the SST turbulence model[C]//Proceedings of the 4th International Symposium on Turbulence,Heat and Mass Transfer.West Redding,2003:625-632.

[23] Ingber L.Adaptive simulated annealing (ASA):Lessons learned[J].Control and Cybernetics,1996,25(1):33-54.

[24] Deb K,Pratap A,Agarwal S,et al.A fast and elitist multiobjective genetic algorithm:NSGA-II[J].IEEE Transactions on Evolutionary Computation,2002,6(2):182-197.

(编辑：吕耀辉)

Aerodynamic optimization of hypersonic lifting body based on FFD method

XIA Chen-chao,SHAO Chun,JIANG Ting-ting,CHEN Wei-fang

(School of Aeronautics and Astronautics,Zhejiang University,Hangzhou 310027,China)

A general method of parametric modeling and grid generation was developed by combining the free form deformation(FFD) and transfinite interpolation techniques.Smoothness of geometry shape and satisfactory quality of grid were achieved from two-and three-dimensional cases.An optimization framework was then constructed by adopting a surrogate model and optimization algorithms.The aerodynamic shape optimization of a hypersonic lifting body was utilized as a demonstration case for the developed framework.The result of single objective optimization using adaptive simulated annealing algorithm shows that the lift-drag ratio has an increase of 1.28% with a constraint on the volume.Besides,multi-objective optimization is carried out by using NSGA-II algorithm,and the Pareto solutions about volume and lift-drag ratio are obtained.One of the typical solutions on the Pareto front shows that the lift-drag ratio has an increase of 2.93% and the volume has an increase of 2.49%.Optimization of the lifting body in this work shows the usability of the free form deformation technique and the good practicality of the optimization framework.

free form deformation;transfinite interpolation;hypersonic;lifting body;design optimization

2014-08-15；

：2014-11-17。

国家重点基础研究发展计划(2014CB340201)。

夏陈超(1989—)，男，博士生，研究方向为飞行器优化设计。E-mail：aeroxia@zju.edu.cn

V441

A

1006-2793(2015)06-0751-06

10.7673/j.issn.1006-2793.2015.06.001