基于全局灵敏度分析方法的气动设计研究

夏露, 杨梅花, 李朗, 张欣

(西北工业大学 航空学院, 陕西 西安 710072)

随着计算机技术的发展,基于智能优化算法与CFD技术相结合的气动优化设计方法已经获得了大量的研究与应用,其优点是具有良好的全局性、适应性和鲁棒性,且使用方便灵活[1]。但对于复杂气动布局而言,设计变量对输出目标的影响差异明显,若不加以区分处理,会对全局优化搜索造成一定的影响,同时,对于机翼、全机等具有大规模设计变量的优化问题,优化设计的复杂程度以及计算量的需求会随着设计变量的增加而急剧上升,对于一般的优化框架在优化效果与优化效率上构成了很大的挑战。针对以上问题,许多学者对全局灵敏度分析方法在气动设计中的应用展开了研究,如:国外学者Sobieski等[2]阐明了灵敏度分析方法在多学科优化设计领域的重要性;Peter等[3]总结了灵敏度分析在气动优化设计中的一些应用方法。国内学者颜力等研究了复变量、自动微分等灵敏度分析方法在多学科设计中的应用[4];姜亚楠等则研究了灵敏度分析方法在飞行器稳健优化优化设计中的应用[5]。由于灵敏度方法在分析和处理输入变量对输出值的影响上具有强大优势,在气动优化设计中也具有明显的优点,本文在前面学者研究的基础上,着力研究灵敏度分析方法在气动优化设计中的应用,以此丰富气动设计手段,提升气动优化效率与优化结果。

本文将全局灵敏度分析方法与气动优化设计相结合,根据设计变量对目标函数的灵敏度信息对设计变量进行分层,然后建立基于全局灵敏度分析方法的分层优化设计系统,对翼型以及机翼进行分层优化设计,并与普通的全参数优化系统对比分析。计算结果表明分层优化系统可以有效减少单次优化时的设计变量数目,减少优化搜索难度,加快优化收敛速度,同时获得较好的优化结果。

1 灵敏度分析方法研究

全局灵敏度分析方法是在可行域内研究输入变量对系统输出的影响,其特点是不仅能分析出单个输入变量对输出的影响效果,同时也能分析输入变量之间的相互作用对输出的影响效果,并且其分析过程具有全局性,能较准确地反应系统特征,因而被广泛应用于各个领域。全局灵敏度分析方法主要分为3类,蒙特卡罗类方法、筛选法和基于方差的方法。蒙特卡罗类方法大多是使用统计学的方法,其计算方法较为简单,但对非线性、非单调问题的灵敏度分析精度并不高。在筛选法中OAT方法计算原理简单,计算效率高,应用范围较为广泛。在方差法中Sobol方法在灵敏度的定量分析方面具有较强的优势,不过存在计算量偏大的问题。下面重点介绍OAT方法和Sobol方法,并利用这2种灵敏度分析方法对相关函数进行灵敏度测试分析,以研究这2种灵敏度分析方法的性能与特点。

1.1 M-OAT方法

筛选法中最简单的一种方法是OAT方法,这种方法的特点是每次只能分析一个因素对输出的影响。OAT方法有很多种,其中Morris于1991年提出的OAT方法使用更为广泛,以下简称为M-OAT方法[6]。

假设设计变量x=[x1,x2,…,xn],共有n个分量。对于某个样本点,xi对输出值的初步灵敏度值可用下式计算:

di(x)=[y(x1,…,xi-1,xi+Δ,xi+1,…,xn)-y(x)]/Δ

(1)

式中，Δ为设计变量的变化量。先将设计变量x的每个分量xi的取值范围归一化与离散化,每个分量从{0,1/(k-1),2/(k-2),…,1}中取值。其中,k为参数的取样点个数。经过以上划分,n维的取值空间Ω就拥有了n×k个取样点,xi便可从这些取样点中随机取样。取样矩阵方程如下:

(2)

矩阵方程中每相邻2行只有一个参数的取值不同, 且其变化量为Δ,把这2行作为模型的参数输入, 分别获得模型的输出结果yi和yi+1然后利用公式di(x)=(yi+1-yi)/Δ就能求出这个样本点中xi的灵敏度值。进行n+1运算就可以得到x1,x2,…,xn的灵敏度值,重复上述的方式计算r次,求取平均值即可得到xi的灵敏度信息Si,即:

(3)

1.2 Sobol方法

Sobol方法是Sobieski在1990年提出的一种基于Monte Carlo积分灵敏度分析方法[7],是一种基于方差分解(analysis of variance,ANOVA)的方法。设所研究的计算模型的输入输出关系是由函数Y=f(X)表示,其中X=(X1,X2,…,Xn)表示n维输入变量,Y表示一维输出变量,设模型Y=f(X)可以分解为如下形式:

(4)

式中:

f0=E(Y)

fi=E(Y/Xi)-f0

fij=E(Y/Xi,Xj)-fi-fj-f0

模型函数Y=f(X)可以分解为维数递增的2n个子函数项。若输入变量之间是相互独立的,则(4)式右端各项相互正交(即协方差为0),因此,可对左右两端求方差得到:

(5)

令Si=Vi/V(Y),则:

(6)

同理:

(7)

以此类推,并对(5)两端同除以V(Y),得到如下公式:

(8)

式中Si表示单个变量Xi对输出方差的贡献效果,也被称为主效应,是Sobol法中的一阶灵敏度指标;Sij表示变量Xi与Xj的交叉影响效果,即二阶灵敏度指标;同理对于S12…m(其中m≤n)表示的是一组输入变量X(X1,X2,…Xs)对输出的交叉影响。

2 算例测试

由于M-OAT和Sobol方法都依赖于样本数据,不同的样本数量会对结果产生重要影响,因此在具体计算过程中,逐渐增加样本数量,根据不同样本数量下的计算结果来判断灵敏度分析方法的分析效率和分析精度。算例中将会用到2个测试函数,如表1所示。2个测试函数中n均代表矢量x的维数,xi为[0,1]之间均匀分布的输入变量,参数ai影响着xi的灵敏度大小,且ai越小,灵敏度值则越大。

表1 测试函数

2.1 函数测试1

选择G1函数作为测试函数,算例中,设n=6,Si为xi对函数的灵敏度值,a=[8,2,0.5,1,4,16], 研究6个输入变量对输出的影响。由参数a的设置可知xi的灵敏度值大小排序理论上应为:

S3>S4>S2>S5>S1>S6

分别使用M-OAT方法和Sobol方法对函数进行设计变量的灵敏度计算,以检验各灵敏度分析方法性能。

图1是M-OAT方法的计算结果显示图,样本数以28为基础,每次增加14个样本。由图1可知,大约在56个样本的时候,MOAT方法就已获得不错的结果,且随着样本数的增加,灵敏度信息变化并不大。由此可知,MOAT方法对样本数量的要求较低,且收敛速度较快。

图2是Sobol方法的计算结果显示图,样本数以32为基础,每次增加16个样本。由图2可知,大约在120个样本的时候,Sobol方法才获得较不错的结果,之后随着样本数的增加,灵敏度信息还有一定的波动。由此可知,Sobol方法对样本数的依赖性较强,当样本数量很少时,计算结果有一定的误差,收敛速度较慢。

图1 M-OAT对G1函数计算灵敏度情况 图2 Sobol方法对G1函数计算灵敏度情况

2.2 函数测试2

本算例选择的测试函数为G2函数。xi为[0,1]之间均匀分布的输入变量。本算例中设n=10,a=[16,8,4,2,1.2,0.7,0.5,0.25,0.1,0],即研究10个设计变量对输出的影响。由参数a的设置可知,输入变量的灵敏度信息排序理论上应为:

S1

图3是M-OAT方法的计算结果显示图,样本数以330个为基础,每次增加330个样本。由图可知,1 650个样本的时候就能获得很好的结果。通过分析1 650个样本和1 980个样本的数据可知,对于同一输入变量,其误差大约在5%左右。说明即使大幅增加样本数量,其分析精度的提升也不明显。

图4是Sobol方法灵敏度分析结果显示图。Sobol方法对样本数的依赖较大,所以本算例中分别取1 980、12 000、16 008、18 000个样本再次进行验证。由图4可知,对于10个输入变量的问题,Sobol方法大约在16 008个样本的时候才获得较好的结果,然而M-OAT在1 650个样本时就获得了不错的结果。通过分析16 008个样本和18 000个样本的数据可知,对于同一输入变量,其误差大约在1%左右,其计算精度达到了相当高的水平。

由以上分析可知,M-OAT方法的特点是当样本数量较少时,其已有相对不错的分析精度,即对样本依赖性较小;而Sobol方法的特点是当样本数量较少时,其计算精度较差,只有当样本数量很大时,其计算精度较高, 即对样本依赖性较大。且根据计算耗时可知,M-OAT方法计算效率较高,而Sobol方法计算效率较低。

图3 M-OAT对G2函数计算灵敏度情况 图4 Sobol方法对G2函数计算灵敏度情况

3 基于全局灵敏度分析的分层优化设计

考虑到气动优化设计较为复杂,计算量偏大,且不需要太高的灵敏度分析精度,本文选用M-OAT方法建立分层优化设计系统。具体方法如下:先利用M-OAT方法对设计变量进行灵敏度分析,再根据设计变量灵敏度大小对其进行分层,然后分别对每一层设计变量进行优化。在对某一层设计变量进行优化时,其余层内的设计变量设为最优定值,当所有层次设计变量优化完毕之后,其最优结果可作为一轮优化的结果。若要充分考虑到高低灵敏度设计变量间的耦合效果,可在一轮优化完成后,更新当前最优解进行下一轮优化,直到优化收敛为止。

3.1 分层优化设计框架

基于以上分析,本文建立的基于全局灵敏度分析的分层优化设计系统把一个复杂的优化问题简化为n个相对简单的子优化问题,其框架如图5所示:

图5 分层优化设计框架

3.2 翼型分层优化设计

本节将使用分层优化设计框架对翼型进行优化设计,并与普通全参数优化框架的结果进行对比,为尽量减小计算复杂程度,分层优化框架和普通全参数框架的算法均使用普通粒子群算法,种群规模为25[8]。对RAE2822翼型进行单个设计点的减阻优化设计。

设计状态:Ma=0.73,α=2.28°,Re=5.51×106

优化目标为最小化翼型阻力系数,约束条件为保持升力系数上下浮动2%,翼型的力矩系数不减小,翼型最大厚度及面积不减小。

约束如下:

式中，CD为阻力系数，CL0为原始升力系数，CL为升力系数，A为翼型面积，A0为原始翼型面积，t为翼型最大厚度，t0为原始翼型最大厚度，CM为力矩系数，CM0为原始翼型力矩系数。

采用CST参数化方法对翼型进行参数化,上下表面各6个控制参数,即设计变量总数为12个。使用雷诺平均N-S方程进行流场计算,湍流模型选k-ωSST模型,采用结构网格,并利用调用脚本文件方式自动生成网格。

3.2.1 灵敏度分析结果

设计变量灵敏度分析结果,如图6所示,其中X1～X12代表翼型的12个设计变量,且以从翼型上表面到下表面,从前缘到后缘的顺序。

图6 各设计变量灵敏度分析情况

由图6可以看出,104个样本之后的灵敏度信息变化不大。根据灵敏度信息值可知第2～4设计变量对阻力系数的灵敏度值较大,而这几个设计变量主要影响着激波的强弱与发展,与实际情况较为符合。综合设计变量数及灵敏度信息,考虑对其进行2层1轮优化,即将灵敏度最大的前5个设计变量分类为第1层,而其余7个设计变量分类为第2层,然后对第1层、第2层设计变量进行一轮优化设计。

3.2.2 翼型优化结果对比

根据上面的灵敏度分析结果,对RAE2822翼型进行分层优化设计。第1层的迭代次数取为19代,CFD调用次数为500次(25+25×19=500),第2层的迭代次数取为20代,CFD调用次数为500次(25×20=500),并与迭代次数为39代的普通全参数优化框架下的结果进行对比,如表2所示。其中Opt1和Opt2分别是第1层和第2层优化结果,OptN是普通全参数优化结果。

表2 普通优化与分层优化气动数据对比

由表2可知,Opt1的阻力系数比OptN的稍大,而Opt2的阻力系数与OptN的结果却非常接近,表明本文建立的分层优化框架可以在单次优化过程中减小设计变量数目的同时获得相对不错的优化结果。图7为分层优化框架与普通优化框架收敛历史对比情况。由图可知,分层优化框架中第1层优化迭代12次接近收敛,第2层优化迭代6次接近收敛,而普通优化需35次才接近收敛。与普通全参数优化设计相比,分层优化框架第1层和第2层优化收敛速度均较快,说明将高低灵敏度设计变量分开优化之后,可以有效减小设计变量之间的干扰,使得优化收敛速度变快。

图8和图9是Opt1、Opt2和OptN的压力系数对比,从压力系数分布对比图可以看出,初始翼型上表面存在一个较强的激波,经过Opt1优化之后,强激波得到削弱,但仍存在一个弱激波;经过Opt2优化之后,弱激波也得到消除;OptN通过39代的优化后,激波同样得到消除,Opt2与OptN相比,压力系数差别不大。

图7 分层优化与普通优化的阻力系数优化历程比较 图8 Opt1与Opt2翼型压力系数对比 图9 Opt2与OptN翼型压力系数对比

图10和图11为原始翼型与各优化翼型对比图。从翼型形状来看,Opt1的变化主要体现在头部有所增大,Opt2在Opt1基础上减小了下表面的头部,前加载有所增强。OptN与Opt2的翼型有所类似,都是在RAE2822基础上增大头部,下表面头部有所减小。

图10 Opt1与Opt2翼型对比 图11 Opt2与OptN翼型对比

由以上分析可知,在单次优化过程中减小了设计变量数目的情况下,分层优化框架依然能获得与普通全参数优化相近的结果。并且与普通全参数优化框架相比,分层优化框架的第1层与第2层优化均能在较快的速度下达到收敛结果,即说明本文建立的分层策略是有效的。由于本算例是针对翼型的优化设计,变量数较少,因此,还需对具有更多设计变量的优化问题进行分析验证。

3.3 机翼分层优化设计

为了检验基于全局灵敏度分析方法的分层优化框架在拥有大量设计变量的气动优化设计中的特性,对机翼进行分层优化设计,并与普通全参数优化进行对比。本算例使用3个翼型剖面来控制整个机翼,取翼根,拐折,翼尖处为控制面,且3个初始剖面选用NACA0012翼型[9]。采用Hicks-Henne参数化方法对翼型进行参数化,上下翼面各6个参数,即1个翼型剖面设计变量为12个,3个控制翼型设计变量总数为36个。并以机翼半展长、前缘后掠角、翼根弦长、拐折处弦长、拐折处展向位置、拐折处翼型扭转角、翼尖弦长和翼尖扭转角8个参数作为设计变量,则设计变量总数为44个。设计状态如下:

Ma=0.78,α=3°,Re=5.89×106

约束状态如下:

式中，troot-max为翼根最大厚度，tbreak-max为拐插处最大厚度，ttip-max为翼尖最大厚度，S为机翼面积。

使用N-S方程进行流场计算,湍流模型选k-ωSST模型,优化算法采用粒子群算法,种群规模为30。

3.3.1 灵敏度分析结果

先使用普通的全参数优化框架对机翼进行优化,粒子群迭代次数取为19。首先对所有设计变量按从翼根到翼尖、从上表面到下表面、从前缘到后缘,再加上机翼半展长等8个参数的顺序依次编号为1～44,使用M-OAT方法对设计变量进行灵敏度分析。分别选取180,360个样本进行灵敏度分析,其灵敏度分析结果如图12所示。

图12 各设计变量灵敏度分析情况

将所有设计变量按照灵敏度值从大到小排序,综合设计变量数目和灵敏度信息,本文考虑对其进行2层2轮优化设计,即将灵敏度最大的前23个设计变量分类为第1层,而其余21个设计变量分类第2层,随后对第1层、第2层设计变量进行2轮优化设计。

3.3.2 机翼优化结果对比与分析

根据上面的灵敏度分层结果,对机翼进行分层优化设计。在第1轮优化设计中,第1层与第2层的迭代次数均取为5代,在第2轮优化设计中,第1层的迭代次数取为5代,第2层的迭代次数取为4代,则CFD总调用次数为600次,并与迭代次数为19次(CFD调用次数为600次)的普通全参数优化框架结果进行对比,其优化结果分别如表3和图13所示。其中1-Opt1和1-Opt2分别是第1轮优化中的第1层和第2层优化结果,而2-Opt1和2-Opt2分别是第2轮优化中的第1层和第2层优化结果,OptN是普通全参数优化结果,Original为原始机翼的数据。

表3 优化前后机翼的部分气动参数比较

图13 分层优化与普通优化框架收敛历程对比

从表3的气动参数数据来看,经过2轮分层优化之后,2-Opt2的阻力系数比OptN的小,表明分层优化框架在单次优化过程中有效减小设计变量数目的同时,由于高低灵敏度设计变量间的信息得到共享,使得分层优化系统能够获得较好的优化结果。

由图13阻力系数收敛历程对比可知,分层优化框架在第1轮优化前期和第2轮优化前期均有较快的收敛速度,说明将高低设计变量分开优化,可以有效减小设计变量间的相互干扰,从而加快优化收敛速度。

图14为优化前后机翼3个翼型控制面归一化形状对比图。

图14 3个翼型控制面归一化形状对比图

通过对翼型和机翼的优化算例分析可知,与普通的全参数优化框架相比,本文建立的分层优化设计系统在单次优化设计过程中能有效减小设计变量数目的情况下,依然能获得较好的优化结果,这大大降低了优化设计的难度。同时由于将具有不同灵敏度的设计变量进行分层优化,可以有效减弱设计变量间的相互干扰问题,从而加快优化收敛速度。

4 结 论

1) 针对气动设计中遇到的具有大规模设计变量的优化设计问题,本文建立了基于全局灵敏度分析方法的分层优化设计系统,与普通全参数优化设计系统相比,分层优化系统可以有效减少单次优化时的设计变量数目,减少了优化搜索难度,加快优化收敛速度,同时获得较好的优化结果。对于拥有大规模设计变量优化设计而言,本文建立的基于全局灵敏度分析方法的分层优化设计方法是一种有效的设计手段。

2) 本文建立的分层优化系统主要针对具有大规模设计变量的优化设计问题,但由于时间、篇幅等因素限制,本文所选算例设计变量数较少,后续工作中考虑大量增加设计变量以检验分层优化系统性能。同时,可以在引入灵敏度分析方法的同时引入代理模型,可以有效提高计算效率。

[1] 夏露, 高正红. 基于PARETO的系统分解法及其在飞行器外形优化设计中的应用[J]. 西北工业大学学报, 2006, 24(1):89-93

Xia Lu, Gao Zhenghong. Improving System Decomposition Method for Application to Multi-Objective Optimization of Aircraft Configuration Design[J]. Journal of Northwestern Polytechnical University, 2006, 24(1):89-93 (in Chinese)

[2] Sobieski Sobieszsczanski, Haftka R. Multidisciplinary Aerospace Design Optimization: Survey of Recent Developments[R]. AIAA-1996-0711

[3] Peter J E V, Dwight R P. Numerical Sensitivity Analysis for Aerodynamic Optimization: a Survey of Approaches[J]. Computers & Fluids, 2010, 39(3):373-391

[4] 颜力,陈小前,王振国. 飞行器MDO中灵敏度计算的自动微分方法[J]. 国防科技大学学报,2006,28(2):13-16

Yan Li, Chen Xiaoqian, Wang Zhenguo. Study of the Automatic Differentiation Method for Calculating Sensitivities in the MDO of Flying Vehicles[J]. Journal of National University of Defense Technology, 2006, 28(2):13-16 (in Chinese)

[5] 姜亚楠,王和平. 基于灵敏度分析的飞行器稳健设计优化方法[J]. 航空计算技术,2009,39(5):48-51

Jiang Ya′nan, Wang Heping. Robust Design Optimization Method to Aircraft Based on Sensitivity Analysis[J]. Aeronautical Computing Technique, 2009, 39(5): 48-51 (in Chinese)

[6] Morris M D. Factorial Sampling Plans for Preliminary Computational Experiments[J]. Technometrics, 1991, 33: 161-174

[7] Sobol′ I M. Global Sensitivity Indices for Nonlinear Mathematical Models and their Monte Carlo Estimates[J]. Mathematics and Computers in Simulation, 2001, 55(1/2/3): 271-280

[8] 李丁, 夏露. 改进的粒子群优化算法在气动设计中的应用[J]. 航空学报, 2012, 33(10):1809-1816

Li Ding, Xia Lu. Application of Improved Particle Swarm Optimization Algorithm to Aerodynamic Design[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(10): 1809-1816 (in Chinese)

[9] 李丁, 夏露. 一种混合粒子群优化算法在翼型设计中的应用[J]. 航空计算技术, 2010, 40(6):66-71

Li Ding, Xia Lu. Application of a Hybrid Particle Swarm Optimization to Airfoil Design[J]. Aeronautical Computing Technique, 2010, 40(6): 66-71 (in Chinese)