一种高效高精度的气动弹性结构优化方法

吕志斌 万志强 / LYU Zhibin WAN Zhiqiang

(北京航空航天大学，北京 100191)

0 引言

飞机是高度综合的现代科学技术的体现，飞机结构设计又是飞机设计的主要阶段[1]。但随着现代飞机尤其是民用客机与运输机追求高性能和低结构重量等方面的要求，飞机结构柔性较大，导致飞机结构变形增大。加之近代工业在复合材料、智能材料等新型材料研究上的突破和应用，使得气动弹性问题在飞机设计中越来越凸显出来[2]。传统的“设计-校核-修改”的结构设计模式越来越难以满足现代飞机对于高性能和低重量的双重要求。

基于此，随着气动弹性学科的发展，以克服气动弹性效应所带来的不利影响的同时尽可能地降低飞机结构重量为目标的气动弹性结构优化技术已经逐渐成为一种主动设计的手段，应用于飞机设计的各个阶段[3]。

然而目前工程上在进行气动弹性分析求解时气动力数据来源仍然是采用低阶面元法，该方法计算效率高但精度低；而CFD(Computational Fluid Dynamics，简称CFD)方法虽然计算精度高，但效率太低，而且难以实现优化过程的迭代计算。当前的气动弹性优化算法则主要集中于单一的优化算法，如可行方向法、遗传算法等。这些算法均有其自身的局限性，如敏度信息不易获得、搜索精度低、收敛速度慢等问题[4]。同时考虑到气动弹性结构优化本身往往具有计算量大、效率低、时间长、设备要求高等特点，难以满足方案设计初期对时间的要求。

针对上述问题，本文提出了采用基于高阶面元法的静气动弹性分析方法来求解约束响应，使用遗传/混合算法作为寻优算法，并引入Kriging(克立格)代理模型来减少计算量和优化时间。最后通过一个成熟的大展弦比机翼为算例，验证本文提出的静气动弹性结构优化方法的可行性与可靠性，从而为飞机方案设计阶段提供一种新的优化方法。

1 基于高阶面元法的静气动弹性分析

1.1 高阶面元法

高阶面元法与低阶面元法一样都是线性方法，都是通过在物面网格上进行布置偶极子或者源来模拟气流流经物体表面的流场，都是对线化后的小扰动位势流方程进行求解[5]：

(1)

(2)

其中，M∞为来流马赫数，ø为扰动位函数。

而与低阶面元法不同的是，低阶面元法在平板气动力网格上的压力点需布置在每个气动网格的四分之一弦线上，高阶面元法则没有类似限制。除可以按照几何外形建立三维的气动力计算模型外，高阶面元法在气动力计算模型的表面网格上布置连续的点源，在结点上布置连续的偶极子，线性分布奇点强度，如图1所示[5]。

图1 网格元素奇点分布示意图

1.2 静气动弹性分析框架

与常用的基于低阶面元法的Nastran软件不同，本文使用的静气动弹性分析方法的气动力数据来源于基于高阶面元法的Panair开源程序，结构计算基于柔度法的Nastran静力分析模块，刚性配平为求解基本的飞行力学平衡方程，载荷导数计算采用差分法。图2所示即为基于高阶面元法的静气动弹性分析方法框架。

图2 基于高阶面元法的静气动弹性分析方法框架

2 代理模型方法

代理模型方法主要包含两方面的内容：一是构建模型的样本点的取样方法，又称实验设计方法；二是进行数据拟合和预测的模型构建方法，又称近似方法[6]。

2.1 拉丁超立方抽样

拉丁超立方取样方法于1979年由McKay提出，该方法最大的优点是不需考虑样本点的维度和取样数目，通过控制取样点的位置可以尽可能避免取样点在小领域内重合，从而实现取样点在设计空间的均匀分布[7]。

假设需要在m维向量空间中抽取n个样本，其具体步骤为：

1)将m维向量空间中的每一维均按策略(一般采用均匀分布的策略)分为互不重叠的n个区间，使得每个区间有相同的被选择概率；

2)从向量空间每一维的每个区间中随机抽取一个点作为该维该空间的值，按照划分的区间共抽取m×n个值；

3)再从每一维中随机抽取2种选择的点组成m维向量，共抽取n个，至此样本点抽取完成。

如图3所示为用Matlab仿真出来的二维向量空间的拉丁超立方取样的样本点分布图，取样点数目为8个。

图3 二维空间的拉丁超立方取样

2.2 Kriging代理模型方法

Kriging模型是一种改进的线型回归模型，其函数表达式中既包含了线性回归部分，又包含了非参数部分[8]：

y=F(x)+z(x)

(3)

式中，F(x)为线性回归部分，具有全局近似的特点；z(x)为非参数随机部分，具有局部偏差的特点。

线性回归部分可以表示为：

F(x)=f1(x)β1+…+fp(x)βp

(4)

按照f(x)形式的不同，可分为常数型、线性型和二次型。本文选择如下所示的二次型，式中β可通过广义最小二乘估计得到：

f1(x)=1

f2(x)=x1,…,fn+1(x)=xn

(5)

非参数部分具有如下的统计学特性：

E[z(x)]=0

Var(z(x)]=σ2

E[z(xi),z(xj)]=σ2G(xj,xi)

(6)

式中，G(θ,xj,xi)为i、j两点的空间相关方程，该空间相关方程也有很多不同种类，对模型构造影响很大。本文取Guass函数作为空间相关函数：

(7)

式中，ndv为样本点的维度，θk是空间相关函数在样本点第j个方向的参数，为待求解的量。

对于一组给定样本点X=[x1…xn]T及其响应值Y=[y1…yn]T，由样本点得到待测点的响应值为：

(8)

式中c是待求权系数向量。根据无偏条件及估计方差最小条件，结合拉格朗日乘子，可求得c：

(9)

具体推导过程见参考文献[8]。

3 优化算法

与其他优化问题一样，气动弹性结构优化问题也可以用如下的公式进行描述，即在设计空间中寻找一组变量使目标函数(一般为结构重量)最小化[9]：

Min.F(v)

(10)

S.T.gj(v)≤0j=1,2,…,ncon

(11)

(12)

其中，式(10)为目标函数，本文为结构重量最小；式(11)为约束条件，包括静气弹约束(如发散速度约束、变形约束等)、动气弹约束(如颤振速度约束等)、强度约束、操稳约束以及气动约束(如阻力约束等)等；式(12)为设计变量的上下限，指定设计空间的值域，又叫边界约束。

3.1 遗传算法

遗传算法是一种借鉴自然选择和遗传进化机理而发明出来的自适应概率搜索算法，具有全局寻优的能力[10]。作为进化算法的一种，遗传算法从初代个体开始，以适应度函数为评价标准，通过遗传算子实现对个体的自然选择和遗传进化，从而不断改良个体直至满足条件。

图4所示为遗传算法流程图[9]，其基本步骤为：

图4 遗传算法流程图

1)确定优化策略，定义优化参数，包括群体规模、交叉概率、变异概率、程序终止条件等遗传参数，选择编码方式、选择方式、交叉方式、变异方式等；

2)随机产生指定规模的初代群体(即第一代父群体)；

3)计算群体中每个个体的目标函数值以及约束响应值，从而计算个体的适应度值，从而对个体进行评估；

4)判断是否满足程序终止条件，如果满足则计算结束，否则继续进行计算；

5)通过基本遗传操作(选择、交叉和变异)产生新一代子群体，用子群体代替父群体形成新一代父群体，重复步骤3和步骤4。

3.2 分形算法

分形算法是一种借鉴树的生长来实现搜索寻优的启发式算法，具有较强的局部搜索能力。图5所示为分形算法流程图，其中选点、分枝、剪枝为分形算法的三个主要操作运算[11]。分形算法的原理很简单，通过选点来确定当前区域的最优解，为分枝和剪枝提供参考依据；通过分枝来产生新的样本点，计算样本点的适应度(为了统一个体评价标准，此处依旧引入遗传算法中的适应度概念代替目标函数值)以判断是否更新当前全局和局部最优解；通过剪枝来缩小可行域范围，使可行域逐渐逼近最优解。在分形优化算法中，分枝提高了算法的精度，剪枝加快了算法的收敛速度。

3.3 遗传/分形算法

遗传/分形混合算法结合了两种算法的优点，遗传算法用于全局搜索以避免优化算法陷入局部最优解，分形算法用于最优个体周围的小范围局部寻优。图6所示为遗传/分形混合算法的流程图，若该流程中不适用分形算法，则图6为遗传算法的流程图。

确定遗传策略和优化参数后，首先使用遗传算法对优化问题进行一次基本遗传操作。选择遗传操作后的最优个体，以该个体为中心确定分形算法的初始可行域以及初始最优解(即为该最优个体适应度)。然后对上述可行域使用分形算法进行寻优，并用分形算法得到的最优个体代替上一次遗传算法得到的最优个体。完成上述操作后重新评估该代群体，判断是否满足算法的终止条件。

图6 遗传/分形混合算法流程图

3.4 静气动弹性结构优化框架

基于前文提到的基于高阶面元法的静气动弹性分析方法、代理模型以及遗传/分形算法，本文建立了一种高效高精度的静气动弹性结构优化框架，如图7所示。在该优化框架中，遗传/分形算法为优化算法(优化时不考虑分形算法即为遗传算法)，代理模型代替基于高阶面元法的静气动弹性分析方法作为个体目标函数和约束响应的求解器。在实施本文提出的优化方法时，首先应该完成代理模型的构建和精度校验，如图8所示。

图7 静气动弹性结构优化框架

图8 代理模型构建与精度校验流程图

4 算例分析

4.1 算例描述

本文所用算例的气动和结构模型如图9和图10所示，机翼半展长16.27 m，翼尖弦长1.62 m。优化时，本文仅对机翼进行优化，而不改变机身的结构重量。

图9 算例气动模型

图10 算例结构模型

表1给出了算例的优化方案，优化状态为飞行高度在11 200 m、飞行马赫数0.8时的定直平飞状态，副翼和方向舵偏角均为0，依靠升降舵和攻角进行配平。

表1 算例的优化方案

为了机翼的蒙皮满足从翼根到翼尖、从后缘到前缘逐渐递减的设计准则，优化时分别使用不同的线性函数来表示蒙皮的厚度变化。图11和图12为初始模型机翼上下表面蒙皮的厚度分布。

图11 机翼上表面蒙皮厚度分度

图12 机翼下表面蒙皮厚度分度

4.2 代理模型构建与精度校验

在优化前，首先需要完成代理模型的构建和精度校核。参照图8的流程首先在机翼的设计变量空间中使用拉丁超立方取样法选取4 000个样本点作为样本输入用于模型构建，使用随机方法选取500个样本点作为测试点用于精度校核，样本点和测试点均采用基于高阶面元法的静气动弹性分析方法进行目标函数和约束响应计算。

(13)

(14)

复相关系数R的值介于0和1之间，其值越大，说明预测值和真实值相关性越大即越接近。一般而言，0.5～1.0即为强相关。相对均方根误差RMSE与响应值无关，值介于0和1之间，其值越小，说明预测精度越高[12]。

图13～图16为校验点的Kriging模型预测结果与高阶面元法的静气动弹性程序计算结果的对比，表 2说明了Kriging模型的预测表现。从图表的结果来看，Kriging模型的预测精度很高，可以完全满足工程应用和科学研究的精度要求。

图13 结构重量预测误差百分比

图14 翼尖位移预测误差百分比

图15 翼尖扭角预测误差百分比

图16 副翼效率预测误差百分比

表2 Kriging模型预测表现

4.3 优化结果对比与分析

本文使用了遗传算法和遗传/分形算法，结合Kriging代理模型对机翼进行了优化，优化共计进行了200代，图17显示了两种算法的目标函数随迭代代数变化的曲线，表3为算例的优化结果。由于算例使用的机翼是一个成熟的机翼，因此可优化的程度并不大。

图17 优化目标函数曲线

从目标函数曲线来看，遗传/分形混合优化算法相比于单一的遗传算法具有更快收敛速度和更好的全局寻优能力，更容易收敛到最优解，且最优解也更接近全局最优解。从算例的优化结果来看，相比于初始模型，结构重量减小，翼尖变形和扭角有一定程度的提高，副翼效率略微减小，但都在约束条件范围内(由于存在约束条件阈值，使得在优化时约束响应能略微超过约束条件)。

表3 优化结果

5 结论

本文建立了一种基于高阶面元法和遗传/分形算法的静气动弹性结构优化算法，并引入Kriging代理模型方法来加快优化时间，降低优化成本。其中，基于高阶面元法的静气动弹性优化方法用于约束响应的求解，遗传/分形算法用于目标函数寻优。最后以一个大展弦比客机机翼为例，通过与遗传算法的对比，说明了本文提出方法的实用性和更强的寻优效率。

值得一提的是，本文在对算例进行优化时，大部分的优化时间花费在了样本数据的计算上，而代理模型的构建和优化计算总耗时不超过2h，也显示了本文的方法在时间上的高效性。