随着计算机硬件水平的提高和智能优化算法的不断完善,基于计算流体力学(computational fluid dynamics,CFD)的气动外形智能优化已成为现代飞机设计的重要手段。为了扩大设计空间和提高设计精度,往往需要大量的优化设计变量。但是设计变量的增加不仅使整个优化设计过程更加复杂,而且会导致优化所需时间大大增加。

近年来,针对减少设计变量的参数化方法研究也逐渐引起关注。Chang等^[1]基于正交函数对NACA系列的传统翼型和超临界翼型进行分析,结果表明采用10个正交函数能够对超临界翼型进行很好的描述;Robinson等^[2]针对超临界翼型采用Gram-Schmidt正交化方法得到一组正交基函数,仅使用前两阶正交基函数描述翼型,从而大大减少了设计变量;Toal等^[3]使用本征正交分解方法(proper orthogonal decomposition,POD)在优化过程中对较优的设计样本进行几何层次上的过滤,重新构建了设计空间,明显缩短了设计周期,但是该方法以损失设计精度为代价,可能会导致优化收敛到局部最优;Ghoman等^[4]基于POD方法开展了外形参数减少的研究,结果表明本方法可以有效地减少设计变量数目并避开典型的缺陷(比如型函数选取的困难、计算效率低等原始参数化方法常见的问题),但是没有对设计空间的变换进行详细讨论。

与此同时,约束的处理对优化设计的效率和精度都有很大影响,目前主要的约束处理方法有——惩罚函数方法、寻找可行解方法、保留解可行性的方法以及混合方法等^[5],但是这些方法在实际应用中存在很多难题,比如传统的罚函数方法中的惩罚系数难以确定、寻找可行解的方法很多仅适用于某一类问题等。目前在气动外形设计中多采用罚函数一类的方法来处理约束,但是如果设计空间中不满足约束的样本比例很大,就会造成计算资源的浪费和设计精度的下降。

综上所述,本文基于POD方法开展了气动外形优化设计过程中的设计空间降维和非法样本点过滤方面的研究。设计算例表明,相对于不作降维和过滤处理的优化设计,本文所提出的方法在保证寻优精度的前提下,有效提高了优化设计效率。

本文中对翼型的描述采用Hicks-Henne参数化方法,上下表面各采用相同数量的Hicks-Henne鼓包^[6],鼓包函数分别如下:

点击浏览公式

式中:b为鼓包高度,p为Hicks-Henne鼓包位置,x为翼型的x轴坐标点位置。本文中参数化翼型的各个鼓包高度变化范围通过翼型所需的变形量来确定。很显然,通过增加设计变量可以扩大优化设计空间,从而得到较好的优化结果,但是会延长设计周期,降低优化效率。

为了提高气动外形优化设计效率,本文采用代理模型代替耗时的N-S方程计算。常用的代理模型主要包括多项式响应面模型、人工神经元网格模型、径向基函数模型及Kriging模型等,其中Kriging模型具有训练样本点处无偏估计、良好的高度非线性近似能力,非常适合作为代理模型使用,目前Kriging模型在工程优化设计领域得到了广泛应用^[7]。

Kriging模型将未知的函数用一个回归函数B(x)和一个均值为零和方差为σ²的高斯随机过程Z(x)组成,因此未知点的函数值为:

点击浏览公式

通过无偏估计和最大似然估计的方法,可以得到:

点击浏览公式

式中:B(x)=f^T,f是已知点优化变量的函数,对于常用的零阶回归函数,f是一维数组,其值为1,为回归参数,通过回归分析可以得到。Z(x)=r^T(x)R^-1(y-f),y是已知样本点的函数值,R是已知样本点处的相关矩阵,r是未知点和已知样本点之间的相关矢量,最常采用的相关函数为Gauss函数,如下所示:

点击浏览公式

式中:θ=(θ₁,θ₂,…,θ_n)^T是空间相关参数矢量,可以通过最大似然估计得到θ。

文献[3, 4]所采用的POD是一种模型降阶方法,基于“快照”的思想,从大量的样本中提取出主要特征(基模态)^[8]。

POD降阶方法的主要思想是寻找一个子空间使所有“快照”｛Y_i,i=1,2,…,N｝(Y_i∈Rⁿ)在该空间中的投影误差最小。若向模态空间的正交投影关系定义为Φ:Rⁿ→Rⁿ,即

点击浏览公式

(5)式等价于寻找各阶基模态｛Φ_j,j=1,2,…,r｝使

点击浏览公式

对于快照集中的所有快照均成立,其中〈·,·〉是两向量的内积。本文选择常用的实数域内欧式空间中的内积和范数,具体表达式如下:

点击浏览公式

点击浏览公式

计算POD基模态共有2种方法:①传统的特征值分解降阶方法,首先构造自相关矩阵并对其仅需特征分解求得各个特征值及其对应的特征向量,通过特征向量和各个“快照”即可求得各阶基模态,而其所对应的特征值大小表征各阶基模态所含“能量”;②奇异值分解(SVD)方法计算各阶基模态。SVD方法不仅计算效率高于特征值分解POD降阶方法,而且在计算高阶模态时比特征值分解方法更为精确^[9],因此本文采用SVD方法,具体过程如下:

假设“快照”集为｛Y_i,i=1,2,…,N｝,首先得到所有快照的均值,如下:

点击浏览公式

通过快照的脉动Y′_i=Y_i-构造矩阵A如下:

点击浏览公式

将A矩阵进行SVD分解得到:

点击浏览公式

式中,U∈R^N×N,∑∈R^N×n,V∈R^n×n,∑仅有对角线元素∑_ii=σ_i,且满足σ₁≥σ₂≥…≥σ_r≥0,其余元素全部为0。这样就可以得到各阶基模态:

点击浏览公式

第i个“快照”可由所有基模态的线性组合得到:

点击浏览公式

通过“能量准则”可以选择基模态数量以近似表达各个“快照”,比如选取p个基模态使前p个基模态的“能量”之和占所有基模态“能量”之和的99.9％以上,其中基模态的“能量”的大小以其所对应的奇异值平方表征,即

点击浏览公式

最终得到各个“快照”的近似解如下:

点击浏览公式

式中,p < r。

采用粒子群优化算法对RAE2822翼型进行气动外形优化设计,设计状态为:Ma=0.729,α=2.31°,Re=6.5×10⁶。为了保证CFD计算的可靠性,计算结果和风洞试验结果的压力分布对比如图 1所示,计算结果对前缘吸力峰值处的小凸起、上表面压力平台区、压力恢复区等吻合较好,对激波的捕捉有所欠缺,基本达到工程所需的精度要求。

图 1 CFD计算结果与风洞试验结果对比

图选项

设计目标和约束如下:

点击浏览公式

式中,C_d为阻力系数,C_l为升力系数,目标函数为升阻比的倒数C_d／C_l。thickness constraints为:前梁位于0.16c处,后梁位于0.6c处,前梁厚度、后梁厚度和最大厚度均不小于RAE2822翼型的相应厚度。对于厚度约束,分别采用惩罚策略、拒绝策略和基于POD方法的设计空间过滤3种措施进行处理。

惩罚措施是对约束进行处理的最一般的方式,是通过对不可行解的惩罚来将约束问题转化为无约束问题。任何对于约束的违反都要在目标函数中添加惩罚项。但是在实际应用中很难确定惩罚函数的形式和惩罚权重,本文中的适应值函数如下:

点击浏览公式

式中:θ₁,θ₂,θ₃为翼型的前梁厚度、后梁厚度和最大厚度,θ_o1,θ_o2,θ_o3为RAE2822翼型相应位置处的厚度。a₁,a₂,a₃为相应的惩罚权重,为了保证惩罚强度,本文采用惩罚权重随优化代数线性递增的方法,a₁,a₂,a₃的初始值均设置为0.2,最后一代的惩罚权重为100。

拒绝策略也称为死亡惩罚,其做法是直接拒绝优化过程中所有不可行解,减小了搜索范围。由于无需设置格外的参数,所以拒绝策略是处理约束最简单的方法。在种群初始化中,为了保证初始种群全部满足约束,首先采用拉丁超立方取样方法(Latin hypercube sampling,LHS)选取大量样本,将其中满足约束的样本作为初始种群;优化过程中生成新的样本之后首先判断其是否满足约束,若满足则加入种群中,否则重新生成新的样本,直到满足约束为止。

为提高设计空间中满足约束的样本比例,对二维翼型从Hicks-Henne参数化空间向POD模态空间的转换进行研究。具体过程如下:

1) 采用LHS方法在Hicks-Henne参数空间取出大量的样本,并判断样本是否满足厚度约束,满足则放入“快照”集中,不满足则舍弃该样本;

2) 采用第2节中介绍的降阶方法从“快照”集中提取出各阶基模态,并求得所有“快照”在各阶基模态上的投影系数;

3) 从投影系数中提取出所有“快照”在各阶基模态上投影系数的上下限,得到快照在模态空间上的投影系数的变化范围。

这样就把优化问题从Hicks-Henne参数化空间转换到模态空间。具体步骤如下:首先采用Hicks-Henne参数化方法和拉丁超立方取样方法^[10],上下表面各选16个鼓包,共有32个参数,得到20 000个翼型,其中满足几何约束的约占24％左右,对这些翼型进行SVD分解得到各阶基模态及其对应的奇异值。图 2是得到的各阶基模态的“能量”的变化情况,可见“能量”主要集中在前若干个基模态中。

图 2 POD各阶模态的“能量”大小

图选项

所选取基向量的“能量”之和超过所有基向量“能量”之和的99.9％,最终确定的基向量数目为16。

多次采用LHS方法在32维的Hicks-Henne参数化空间和16维的POD基模态空间分别获得随机不同的20 000个翼型样本发现,结果如图 3所示。由图可见,Hicks-Henne参数空间中满足厚度约束的样本比例约占24％,而在POD基模态空间中该比例则超过70％。由此可见,从参数化空间转换到基模态空间之后,设计空间有所减小,许多非法的设计空间被剔除,有利于缩短优化耗时。

图 3 不同设计空间中取样结果对比

图选项

为了体现POD空间重构方法在优化设计中的可行性,本文采用32个Hicks-Henne鼓包构成的参数化空间和前16个POD基向量对应的设计空间分别对RAE2822翼型在厚度约束下的升阻比进行优化。

在气动外形优化设计之前,必须首先训练具有高可信度的Kriging代理模型,从而减少优化设计耗时,提高效率。在Hicks-Henne参数化空间和POD基模态空间中分别采用LHS方法得到500个和300个训练样本,以此训练Kriging代理模型。

然后从两个空间中随机抽取50个测试样本,使用已建立的代理模型对测试样本的升力和阻力系数进行预测。在Hicks-Henne参数化空间和POD基模态空间中建立的代理模型对测试样本升力系数预测均方根误差分别为1.25％、1.46％,阻力系数预测误差分别为2.77％和3.77％。其中升力系数的计算结果和预测结果对比如图 4所示。

图 4 计算和预测的结果对比

图选项

在优化过程中适当地对代理模型进行更新以保证其精度和可靠性,本文采用的方法是每代选取3个适应值最高的样本,采用CFD计算得到其气动力系数,并将其加入代理模型样本库,对其重新进行训练,得到新的Kriging代理模型。

为了对比翼型优化设计在Hicks-Henne参数化空间和POD基模态空间中的不同表现,采用粒子群智能优化算法对翼型进行优化,优化代数为50代,每代种群由100个样本构成。每种策略(Hicks-Henne参数化空间+罚函数、Hicks-Henne参数化空间+拒绝策略、POD基模态空间+罚函数、POD基模态空间+拒绝策略)优化5次,表 1为每次优化得到最优构型的适应值对比。

由表 1可见,16个POD基向量系数的优化结果略优于32个Hicks-Henne参数化方法的优化结果,而拒绝策略的结果也略好于罚函数方法。其中每种策略的最优收敛过程对比如图 5所示。

图 5 收敛结果对比

图选项

采用12核Intel(R) Xeon(R) CPU X5650 @2.67GHz图站开展优化设计,对比采用拒绝策略的Hicks-Henne方法和POD方法的优化耗时如表 2所示,表中时间单位为h。优化过程中的优化耗时主要用于Kriging代理模型的更新。更新代理模型过程中需要大量调用BLAS库函数中的向量运算程序,较为耗时。相比于Hicks-Henne方法,采用POD方法节省了约1／3的优化耗时。

最终优化结果的几何和压力分布对比如图 6和图 7,升阻力系数见表 3。4个构型的升阻比相差在2.0以内,前缘吸力峰值压力系数差量在0.1以内,激波均有了明显减弱。其中采用Hicks-Henne方法优化结果在翼型上表面70％至90％弦长位置有明显的凸起,该凸起能够降低逆压梯度,采用罚函数策略的最为明显,该优化翼型上表面基本无激波;采用POD方法得到的2个最优构型均有明显的前加载特征;而4个优化结果在翼型下表面后缘处的后加载强度相差不大。

图 6 优化前后几何外形对比

图选项

图 7 优化前后压力分布对比

图选项

1) 对于厚度约束情况,Hicks-Henne参数化方法得到的设计空间中满足约束的样本约占24.36％,而基于POD方法重构的设计空间中满足约束的样本约占70.26％,大大提高了合法样本的比例;

2) 在基本不降低设计精度的前提下,基于本征正交分解方法的气动外形优化设计空间重构能够有效地减少了气动外形优化设计参数,缩短设计周期,提高寻优效率;

3) 目前最常用的2种惩罚策略——罚函数法和拒绝策略,罚函数方法要反复试验方能得到较好的惩罚因子,从而提高设计精度,但是拒绝策略较为简单。在厚度约束下RAE2822翼型单目标增加升阻比优化中拒绝策略得到的优化结果略优于罚函数法。